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We present a realization of astronomical relativistic reference frames in the Solar System and its 
application to the GRAIL mission. We model the necessary spacetime coordinate transformations 
for light-trip time computations and address some practical aspects of the implementation of the 
resulting model. We develop all the relevant relativistic coordinate transformations that are needed 
to describe the motion of the GRAIL spacecraft and to compute all observable quantities. We take 
into account major relativistic effects contributing to the dual one-way range observable, which 
is derived from one-way signal travel times between the two GRAIL spacecraft. We develop a 
^s^j , general relativistic model for this fundamental observable of GRAIL, accurate to 1 /im. We develop 

and present a relativistic model for another key observable of this experiment, the dual one-way 
range-rate, accurate to 1 /im/s. The presented formulation justifies the basic assumptions behind 

Qthe design of the GRAIL mission. It may also be used to further improve the already impressive 
results of this lunar gravity recovery experiment after the mission is complete. Finally, we present 
qq ' transformation rules for frequencies and gravitational potentials and their application to GRAIL. 
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I. INTRODUCTION 



Several past, present and planned space missions utilize a pair of spacecraft orbiting a celestial body in a tight 
5jT)- formation. Continuous high-precision range and range-rate measurements between the spacecraft yield detailed in- 
formation about the gravity field of the target body. Missions of this type include the Gravity Recovery and Climate 
Experiment (GRACE) mission [1] in orbit around the Earth; the Gravity Recovery and Interior Laboratory (GRAIL) 
mission, which comprises two spacecraft in orbit around the Moon [2|-|5(; an d planned missions such as a GRACE 
£Ni . Follow-on mission or a proposal for a GRAIL-like mission in orbit around Mars, 
(f) ■ Of these, the mission of particular current interest is GRAIL, as the two GRAIL spacecraft are presently (2012) 
! orbiting the Moon. In this paper, we therefore focus on the GRAIL mission and its science observables. However, the 
■ lessons learned are also applicable to other, similar experiments, 
^vq | To reach its science objectives, the GRAIL mission relies on precision navigation of both spacecraft and accurate 
t— I • range measurements between the two lunar orbiters performed with their on-board Ka-band ranging (KBR) system. 
The instantaneous one-way range measurements performed at each spacecraft are time-tagged and processed on 
the ground to form dual one-way range (DOWR) measurements [6|. The mission relies on precision timing of all 
critical events (using the on-board ultra-stable oscillator, or USO) related to the transmission and reception of various 
microwave signals used on GRAIL for formation tracking and navigation. The resulting time series of highly accurate 
radio-metric data will allow for a major increase in accuracy when studying the gravity field of the Moon. The 
differential nature of the science measurements allows for the removal of a number of measurement errors introduced 
in the process. In particular, the approach compensates for errors due to long-term instabilities of the on-board USOs. 
This allows for an improvement in accuracy by about two orders of magnitude when compared to other techniques. 
In fact, the anticipated accuracies are of the order of 1 (im in range and 1 /im/s in range rate. 

It was recognized early on during the mission development that due to the expected high accuracy of ranging 
data on GRAIL, models of its observables must be formulated within the framework of Einstein's general theory of 
relativity. In fact, a naive application of the observable models developed for the GRACE mission yjj may have led 
to significant model discrepancy (as emphasized in Ref. [fjj]), as these models do not take into account relativistic 
contributions that are critical for GRAIL. The ultimate observable model for GRAIL must correctly describe all the 
timing events occurring during the science operations of the mission, including both the navigation observables (S- 
and X-band, ~ 2 GHz and ~ 8 GHz correspondingly) and inter-spacecraft tracking (Ka-band, ~ 32 GHz) data. 

The model must represent the different times at which the events are computed, involving the time of transmission 
of the Ka-band signal at one of the spacecraft, say GRAIL-A, at tAo, and the reception of this signal by its twin, 
GRAIL-B, at time £b- In addition, the model must include a description of the process of transmitting S-band and 
X-band navigation signals from either spacecraft and reception of this signal at a Deep Space Network (DSN) tracking 
station at time tc- 
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FIG. 1: Representative geometry (not to scale) of the vectors involved in the computation of the GRAIL observables. "SSB" is 
the Solar System barycenter, "E" is the center of the Earth, "M" is the center of Moon, "EMB" is the Earth-Moon Barycenter. 
"A" and "B" are the positions of the GRAIL-A and GRAIL-B spacecraft, respectively, and "C" is the position of the DSN 
tracking antenna on the surface of the Earth. 



We model the range Rab = |Rab| between the two spacecraft A and B (see Fig. Q] for geometry and notations) as: 



where xem is the vector connecting the Solar System barycenter (SSB) with the Earth-Moon barycenter (EMB), xm 
is the vector from the EMB to the Moon's (M) center of mass, xa and xb are vectors connecting the SSB with the 
positions of the two GRAIL orbiters and vectors y A and y B connect the Moon's center of mass with the orbiters. 

For navigation purposes, both orbiters maintain communication links with a ground-based DSN antenna. The range 
Rac — |Rac| between a GRAIL spacecraft (GRAIL-A, for instance) and a ground-based antenna can be modeled 
as: 



where xe is the vector from the EMB to the geocenter (E), xc is the vector connecting the SSB with the ground 
antenna whereas the vector y c determines the geocentric position of the ground antenna's reference point. 

For actual computations, we use several different reference systems 1 . The Solar System Barycentric Coordinate 
Reference System (BCRS) has its origin at the SSB. The origin of the Geocentric Coordinate Reference System 
(GCRS) is the Earth's center of mass. Positions of DSN ground stations are given with respect to another terrestrial 
coordinate system, the Topocentric Coordinate Reference System (TCRS; see also Ref. [8]). We also consider the 
Lunicentric Coordinate Reference System (LCRS; for additional discussion, see Ref. @), the origin of which is fixed 
at the Moon's center of mass. Finally, we attach to each spacecraft its Satellite Coordinate Reference System (SCRS; 
for a similar approach aimed to construct a reference frame for the GAIA project, see Ref. [13). (We discuss these 
reference frames and their relationships in depth in Sec. |TTJ) 

Equations ((T|) and offer a good starting point to develop an appropriate relativistic formulation for the experi- 
ment. The six vectors involved in Eqs. dTJ— ([2]) can be expressed in terms of their respective points of origin: e.g., xem 
would be expressed in the BCRS, y A and y B in the LCRS, Rab and Rac in the SCRS of GRAIL-A, etc. Each of 
these coordinate systems has a corresponding time coordinate. To compute the vector sums and differences, all vectors 
involved must be converted to a common relativistic space-time reference system. Although in general relativity one 
can introduce any reference frame to describe the experiment, the best practical choice is offered by some realization 



1 Following Refs. [tI |27|. we use the term "reference system" to describe a purely mathematical construction, while a "reference frame' 
is a physical realization of such. 



-Rab = | Rab I = |*b - x A | = |(x EM + x M + y B ) - (*em + x M + y A )| 5 



(1) 



Rac = | Rac I = |x c - x A | = |(x EM + x E + y c ) - (x E m + x M + y A )|, 



(2) 
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of the BCRS. We will use a realization of the BCRS that is called the SSB reference frame. The coordinate time 
associated with the BCRS is TCB (Barycentric Coordinate Time). For practical applications, it is often preferable 
to use another time scale, the TDB (Barycentric Dynamical Time). Currently published planetary ephemerides are 
provided using TDB. TDB and TCB differ only by a linear scaling. The advantage of using TDB is that the difference 
between it and terrestrial timescales (e.g., TT, defined in Sec. Ill D| ) is as small as possible and periodic. The choice of 
the TDB as the SSB time coordinate is realized by the appropriate linear scaling of space coordinates and planetary 
masses (see [Il|-[l3[ for review). 

The vectors xe,xm, and xem are readily available in the SSB reference frame, obtained by numerical integration 
and from Solar System ephemerides [I3|. The vectors y A , y B and y c have to be transformed to the SSB frame from 
geocentric and lunicentric reference systems, respectively. Clearly, the required conversion between reference systems 
also involves conversion of the relativistic time coordinate. The equations of motion of the Moon and Earth, including 
all the relativistic effects at an accuracy even exceeding that of the GRAIL experiment, have already been discussed 
elsewhere here we concentrate on the computation of observables. 

In this paper, we focus on the formulation of a relativistic model for computing the observables of the GRAIL 
mission, with results that are applicable to other past and planned missions with similar observables. We address 
some practical aspects of the implementation of these computations. In Sec. [IT] we discuss all relevant relativistic 
four-dimensional reference systems and the transformations that are required to make the vector sums in Eqs. ([1]) and 
© computable. In Sec. [En] we discuss the process of forming the inter-satellite Ka-band range (KBR) observables 
of GRAIL and derive a model for the dual one-way range (DOWR) observable. We also develop a relativistic model 
for another fundamental observable on GRAIL: the dual one-way range-rate (DOWRR). We conclude with a set of 
recommendations and an outlook in Sec. IIVI 

In order to keep the main body of the paper focused, we chose to present some calculational details in the form 
of appendices. In Appendix [X] we present some important derivations: In Appendix IA II we derive the solution for 
the post-Minkowskian space-time in general relativity, in Appendix I A 21 we derive analytic expressions to describe the 
phase of an electromagnetic signal in gravitational field, and in Appendix I A 31 we discuss the coordinate gravitational 
time delay. In Appendix [B] contains a discussion on the evaluation of the integral that is needed to assess the full 
accuracy of the DOWR observable. Finally, in Appendix [C] we briefly address the transfer of a precision frequency 
reference between the spacecraft and a ground station. 

The notational conventions used in this paper are as follows. Latin indices from the beginning of the alphabet, 
a, 6, c, are used to denote Solar System bodies. Latin indices from the second half of the alphabet (m, n, ...) are 
space-time indices that run from to 3. Greek indices a, (3, ... are spatial indices that run from 1 to 3. In case of 
repeated indices in products, the Einstein summation rule applies: e.g., a m b m — J2m=o a ™b m - Bold letters denote 
spatial (three-dimensional) vectors: e.g., a = (ai,a 2 ,a3),b = (61,62,63). The dot is used to indicate the Euclidean 
inner product of spatial vectors: e.g., (a-b) = a\b\ +0262 + 0363. Latin indices are raised and lowered using the metric 
g mn . The Minkowski (flat) space-time metric is given by j mn — diag(l, —1, —1, —1), so that 7 /iI /a M 6 I/ = —(a • b). We 
use powers of the inverse of the speed of light, c , and the gravitational constant, G as bookkeeping devices for 
order terms: in the low-velocity (v <C c), weak-field (GM/r <^ c 2 ) approximation, a quantity of 0(c~ 2 ) ~ 0(G), 
for instance, has a magnitude comparable to v 2 /c 2 or GM/c 2 r. The notation 0(a k ,b e ) is used to indicate that the 
preceding expression is free of terms containing powers of a greater than or equal to k, and powers of 6 greater than 
or equal to i. 

II. SPACE-TIME REFERENCE FRAMES AND TRANSFORMATIONS 

The theory of general relativity is generally covariant. In the Riemannian geometry that underlies the theory, 
coordinate charts are merely labels. One may choose an arbitrary coordinate system to describe the results of 
a particular experiment. Space-time coordinates have no direct physical meaning and it is essential to construct 
physical observables as coordinate-independent quantities. 

On the other hand, some of the available coordinate systems have important practical advantages. These systems 
are usually associated with a particular celestial body, ground-based facility or spacecraft, thereby yielding a material 
realization of a reference system to be used to describe the results of precision experiments. In order to interpret the 
results of observations or experiments, one picks a specific coordinate system that is chosen for the sake of convenience 
and calculational expediency, formulates a coordinate picture of the measurement procedure, and then derives the 
observable. It is also known that an ill-defined reference frame may lead to the appearance of non-physical terms 
that may significantly complicate the interpretation of the data. Therefore, in practical problems involving relativistic 
reference frames, choosing the right coordinate system with clearly understood properties is of paramount importance, 
even as we recognize that in principle, all (non-degenerate) coordinate systems are created equal 0j- 

In a recent study [l5j . we presented a new approach to investigate the dynamics of an isolated, gravitationally 
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bound astronomical TV-body system in the weak field, slow-motion approximation of the general theory of relativity. 
Celestial bodies are described using an arbitrary energy-momentum tensor and assumed to possess any number 
of internal multipole moments. Using the harmonic gauge conditions together with a requirement for preserving 
conservation laws, we were able to construct the relativistic proper reference frame associated with a particular body. 
We also were able to determine explicitly all the terms of the resulting coordinate transformations and their inverses. 
In this paper we rely on the results obtained in Refs. [lU HH and develop a set of coordinate reference frames for 
GRAIL. 

To reach its scientific objectives, in addition to the BCRS, GRAIL will have to utilize a set of several fundamental 
coordinate reference frames. These include terrestrial reference systems, namely the GCRS and the TCRS, and 
lunar reference systems, the LCRS and SCRS. In Ref. [la ], we presented the detailed structure of the representations 
of the metric tensor corresponding to the various reference frames involved, the rules for transforming relativistic 
gravitational potentials, the coordinate transformations between the frames and the resulting relativistic equations of 
motion. The accuracy that is achievable by these calculations is sufficient to accommodate modern-day experiments 
in the Solar System and exceeds that needed for GRAIL. Here, we present the essential part of these transformations 
between various coordinate systems involved and dealing with transformations of relativistic time scales and position 
vectors, at the level of accuracy required by GRAIL. 



A. Barycentric Coordinate Reference Frame (BCRS) 

The Barycentric Celestial Reference System (BCRS) is defined with coordinates {x m } = (ct,x — x a ), where t is 
TCB. The BCRS is a particular implementation of a barycentric reference system in the Solar System. The metric 
tensor g mn {x) of the BCRS satisfies the harmonic gauge condition. It can be written [15[ as 

.9oo = 1 - -^w + -^w 2 + 0(c~ 6 ), g 0a = -^j a \w x + 0(c~ 5 ), g aj3 = j afj + Jap-^w + C(c~ 4 ), (3) 

where w and w x are harmonic gauge potentials that can be presented, at the level of accuracy suitable for the purposes 
of the GRAIL mission (i.e., neglecting higher order mass- and current-multipole moments), in the form [7L [TM Il6j|: 

w = + * { 2v 2_J2^- Un b -v b f-Hn ■*»)}) +0(c~% (4) 

b b C c^b Tch 

w = £^v b + 0(c- 2 ), (5) 

b 

where = x — z b , = \rj,\, and rib = ffc/^h, with Zb being the barycentric position of body b, and we use r a b = r^ — r a 
to denote the vector separating two bodies a and b. Also, the overdot denotes ordinary differentiation with respect 
to t, Vb = Zb (vb = Wb\) and &b = Zb (&b = \ a b\) are the barycentric velocity and acceleration of body b, and Mb 
is its rest mass. Lastly, the summation in (@|-([5]) is being performed over all the bodies b — 1,2..., N in the Solar 
System. The metric tensor ^ and the gravitational potentials ((U)-© have sufficient accuracy for modern precision 
experiments in the Solar System 0, fl5| . 

From Fig. [T] we can read off the barycentric positions of the Earth, ze, and the Moon, zm- ze = xem + xe 
and zm = xem + xm, respectively. Both of these vectors, and corresponding velocities ve = ze and zm = zm can 
be computed in the first pos t- Newtonian approximation using the Einstein-Infeld-Hoffmann (EIH) equations in the 
coordinates of the BCRS [H H2HII : 



E^{ 1 +?(- 4 E^-E^+^+^- 4 (*-^)-i(^^) a + i('«*- if »))} + 

b^a ab c^a ac c^b bc 



where r a b — \^ a b\ and n a b = r a b/r a b- When describing the motion of spacecraft in the Solar System, the models 
also include forces of attraction between the zonal harmonics of the bodies of interest and forces from asteroids and 
planetary satellites (see details in Ref. (2lT|). 

To determine the orbits of planets and the spacecraft, one must also describe the propagation of electromagnetic 
signals between any two points in space. The light-time equation corresponding to the metric tensor ([3]) and written 
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to the accuracy sufficient for GRAIL has the form (see also Ref. [H, [Hj]): 

. , |r 2 -ri| ^ GM b \ r\ + r b 2 + r b 12 ] 5 

6 6 L L r l + r 2 ~ r 12 J 

where ii refers to the signal transmission time and t 2 refers to the reception time, while i"i i2 are the barycentric 
positions of the transmitter and receiver. Also, r\ 2 are the distances of the transmitter and receiver from the body b 
and r\ 2 is their spatial separation [l^, [2(| ■ The logarithmic contribution in (J7J is the Shapiro gravitational time delay 
that, in the case of GRAIL, is mostly due to the Moon, the Earth, and the Sun. (Note that the 0(c~ 5 ) terms are 
beyond GRAIL's sensitivity; see the analysis in Sec. IIII 51 ) 

The general relativistic equations of motion ([6]) and light-time equation (J7]) are used to produce numerical codes for 
the purposes of constructing Solar System ephemerides, spacecraft navigation [l8l . [2l| and analysis of gravitational 
experiments in the Solar System [20|, l24[. GRAIL also relies on these equations to compute its range and range- 
rate observables between the two spacecraft in lunar orbit. The numerical algorithm developed for this purpose [3| 
iteratively solves the light-time equation ([7]) in the SSB frame in terms of the instantaneous distance between the 
two spacecraft. Our objective is to develop an explicit analytical model for all the quantities involved in these high- 
precision computations. For this purpose, we need a clearly defined set of astronomical reference frames, which we 
discuss next. 



B. Relativistic coordinate transformations between various reference frames 



To describe the dynamics of an iV-body system (such as the Solar System) in general relativity, one may choose to 
introduce N + 1 reference frames, each with its own coordinate chart. We need one global coordinate chart defined for 
the inertial reference frame that covers the entire system under consideration (e.g., BCRS). In the immediate vicinity 
of each of the N bodies in the system we can also introduce a set of local coordinates defined in the frame associated 
with this body (body-centric system). In the remainder of this paper, we use {cc m } to represent the coordinates of 
the global inertial frame and {y™} to be the local coordinates of the accelerated proper reference frame of body a. 

In Ref. [HI, we showed that the transformations between the harmonic coordinates of the BCRS {x m } and non- 
rotating body-centric reference systems {y™} (such as the GCRS or LCRS) may be written in the following form: 



x° = y° a + c- 2 < 
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2"ii TU oxt 






2 V a( V a ' y a 


) - y a u2* - 



dy' a }+O(c- 4 ), 



* = y a +z a + c 2 {iv a (v a -yJ-y Q [/ c a xt + [u; a x y Q ] + ±a Q ^ - y a (y a • a a )} + 0(c 4 ), 



(8) 
(9) 



where z a is the vector that connects the origin of the {x m } reference system with the origin of the {y™} reference 
system. Note that the accuracy of timing for GRAIL is limited by the performance of the on-board USO, which have 
an error of 0(1O -13 ) for 10 3 s of integration time [6j. Therefore, the c -4 terms in Eq. (|8]), which are at most of order 
~ u 4 /c 4 ~ 10~ 16 , arc negligible for GRAIL, even in the absolute sense. The differential nature of the observables 
on GRAIL further reduces the sensitivity of the mission to such small terms in the transformations. For a complete 
post-Newtonian form of these transformations, including the terms c~ 4 and their explicit derivation, consult Ref. [l5j . 
The inverses of the transformations ©-© can be written as 



y° a 



y« 



'{< 
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±v 2 - 

2 a 





dx' 



! {^v Q (v Q ■ r„) + r a C/ c a xt + K x rj + r a (r Q • a a ) - ia^} + 0{c 



(10) 



11) 



where r a = x— z a . The quantity J7" xt in Eqs. ©-© and (jTU]) - (fTTj) is the Newtonian gravitational potential (including, 
if necessary, multipole corrections) due to all bodies in the Solar System other than body a, at the location of body 
a. Furthermore, a a is the Newtonian acceleration of body a due to the combined gravity of all other bodies. Later in 
this section, we will present the expressions for J7° xt and a a for each of the chosen reference frames. 

Finally, uj a in Eqs. ([9]) and (fTTj) is the vector associated with the relativistic precession given as cj" = ^e"to^ v , 
with e" being the fully antisymmetric Levi-Civita symbol, normalized as e3,3 = 1j an( i the matrix uj^ having the 
form [l5i: 



Un' = - 



E 



GMh 



' ba 



{n£,(f«£-2t£) 



n ba \ 2 U a 



2<)}+0( C - 2 ). 



(12) 
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The expression for the relativistic precession matrix is given here only for the sake of completeness. Because of their 
small magnitude (~ 10~ 15 m), these terms will not affect the GRAIL measurements (see discussion in Sec. Ill F() . 

In the rest of this section, we discuss four fundamental body-centric reference frames that are useful for collection 
and interpretation of GRAIL data. 



C. Coordinate systems used in the vicinity of the Earth 



In the vicinity of the Earth, two standard coordinate systems are utilized: the Geocentric Coordinate Reference 
System (GCRS), centered at the Earth's center of mass is used to track orbits in the vicinity of the Earth. The 
positions of objects on the surface of the Earth, such as DSN ground stations, are usually given in the Topocentric 
Coordinate Reference System (TCRS). 



1. Geocentric Coordinate Reference System (GCRS) 

When constructing a body-centric coordinate reference frame for a body a at the level of accuracy anticipated for 
GRAIL, it is sufficient to consider only monopole contributions to the external potential of all the bodies in the 
Solar System (for the Earth it is mostly the Sun and the Moon) excluding body a itself 15]. Thus, for the GCRS, 
the Newtonian potential of the external bodies (i.e., excluding the Earth) Ug xt and the corresponding acceleration a a 
that are present in the coordinate transformations Eqs. have the form (l5j : 

= Y. — + 0{c-% a E = - VE& = - J2 GM b T -^ + 0(c~% (13) 

b#E b#E r bE 

where summation is performed over all the bodies excluding the Earth (symbolically, b ^ E), the vector that connects 
body b with the Earth's center of mass is represented by r bE = xe — x b and the contributions of the higher multipole 
moments of mass distribution within the bodies are neglected due to their smallness. 

The transformations given by Eqs. ([H])-®, together with the potential Ug xt and the acceleration a E given by 
Eq. (|13p . determine the metric tensor g^ n of the non- rotating GCRS [HI . We denote the coordinates of this reference 
frame as {y^} = (j/e^e) an d present the metric tensor gf nn in the following form: 

2 2 4 2 

9oo = 1 - ^ W [E] + -^ W [E] + 0(c -6 ), 9o a = -la\-^w^ + 0{c~ 5 ), g aj3 = 7a/3 + 7a/3^u> E + 0(c~ 4 ), (14) 

where u> E and are the scalar and vector harmonic potentials that are given by 

w E = UE + u^ + Oic- 4 ), (15) 

w E = 4[y E xS E ]+0( c - 2 ), (16) 
z Ve 

where S E in Eq. (|16[) is the Earth's angular momentum. The scalar potential we is formed as a linear superposition 
of the gravitational potential Ue of the isolated Earth and the tidal potential it E ldal produced by all the Solar System 
bodies (excluding the Earth itself, b ^ E) evaluated at the origin of the GCRF. The Earth's gravitational potential 
Ue at a location defined by spherical coordinates (j/e, <fr, 0) is given by 

Ue = ^^{l + ^Y.{—)" P ^ cose ^ C fk^ k( t> + Sf k ^^ (17) 

^ E 1=2 fe=0 ^ E 

where i?oE being the Earth's radius, Pgk are the Legendre polynomials, while Cf k and Sf k are spherical harmonic 
coefficients that characterize the Earth. At the level of sensitivity of GRAIL, only the lowest order spherical harmonic 
coefficients need to be accounted for, and time-dependent contributions due to the elasticity of the Earth can be 
ignored. Insofar as the tidal potential u E ldal is concerned, for GRAIL it is sufficient to keep only its Newtonian 
contribution (primarily due to the Moon and the Sun) which can be given as usual: 

M tMai = £ f Ub(rm + yE) _ Ub(rbE) _ yE . V [/ b (r hE )) ~ £ ^(3(n bE ■ y E ) 2 - y|) + e% E ,^ 2 ), (18) 

b^E b^E r&E 
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where Ub is the Newtonian gravitational potential of body b, r bE is the vector connecting the center of mass of body b 
with that of the Earth, and V denotes the divergence with respect to y E . Note that in Eq. (fTgf we omitted relativistic 
tidal contributions of 0(c~ 2 ) that are produced by the external gravitational potentials. These are of the order of 
10 -16 compared to Ue and, thus, completely negligible for GRAIL. In addition, we present only the largest term in 
the tidal potential of the order of ~ y E ; however, using the explicit form of the tidal potential Eq. (p~8|) . one can easily 
evaluate this expression to any order needed for a particular problem. 

The proper time at the origin of the GCRS is called the Geocentric Coordinate Time (TCG), denoted here as £tcg ■ 
It relates to the barycentric time TCB t as 



dtrCG_ 
dt 



1 

1 - — 



E 

b=£E 



GMh 



TbE 



C(c" 4 ) w 1 - 1.48 x 10" 



The Earth's barycentric velocity ve and position ze can be computed from Eq. ©. 



2. Topocentric Coordinate Reference System (TCRS): proper and coordinate times 



(19) 



.rnm 



To obtain the metric of the topocentric coordinate reference system, the TCRS, one can transform the metric 
of the GCRS using coordinate transformations given by Eqs. (H])-©, where the "external" potential U^p xt is the 
gravitational potential w E given by Eq. (fTS)) and evaluated at the surface of the Earth: 

(20) 
(21) 



Ug t = w E (y c ) = U E (y c ) + J2 (u b (r bE + Yc) " U b (r bE ) - y c • VU b (r bE )) + 0( C - 2 ), 



b^E 



a c = -VUg t = -VU E (y c ) - £ (VU b (r bE + y c ) - VC/ fc (r bE )J + 0(cT 2 ), 

b^E 



where y c is the position vector of the DSN station in the GCRS. Note that L^E(yc) must be treated as the potential 
of an extended body and include a multipolar expansion with sufficient accuracy, taking into account time-dependent 
terms due to tidal effects on the elastic Earth. 

The proper time re, kept by a clock located at the GCRS coordinate position yc(t), and moving with the coordinate 
velocity vqo = dy c / dt^co = [^e x Ycli where f2 E is the angular rotational velocity of the Earth at C, is determined 
by 



drg 



i-4 



2 V C0 



^E(y c ) + E yr^ (3(n hE ■ y c f - y c ) + (a E ■ y c ) 



3 „-4\ 



(22) 



where a E is the Earth's acceleration in the BCRS, Eq. (|T3l) . and n^E is a unit spatial vector in the body-Earth 
direction, i.e., n hE = rfcE/l r 6E|, where r bE is the vector connecting body b with the Earth. The term within the square 
brackets in Eq. (f22|) is the sum of Newtonian tides due to the Sun, the Moon, and other bodies at the clock location 
y c . These terms are small for Earth stations (of order 2 x 10~ 17 ) and are negligible for GRAIL. The last term is due 
to non-inertiality of the GCRS and accounts for the Earth's finite size. This term is evaluated to be of the order of 
4.2 x 10~ 13 , which is about 10~ 3 smaller compared to the gravity potential on the surface of the Earth and, thus, it 
can be omitted. 

Therefore, at the accuracy required for GRAIL, it is sufficient to keep only the first two terms in Eq. (|22p when 
defining the relationship between the proper time tq and the coordinate time £tcg : 



G^TCG 



1 



l r 



Ko + ^(y c ) +0(c- 4 ) 



(23) 



At the level of accuracy required for GRAIL, it is important to account in Eq. (l23t for the oblateness (non-sphericity) 
of the Earth's Newtonian potential, which is given in the form of Eq. (|17[) . In fact, when we model the Earth's gravity 
potential, we need to take into account quadrupole and higher moments, time-dependent terms due to tides as well 
as the tidal displacement of the DSN station. For example, for a clock situated on the surface of the Earth, the 
relativistic correction term appearing in Eq. (|23|) is given at the needed precision by 



^ + U E (y c ) =W - £ C gdh, 



(24) 
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where Wq = 6.2636856 x 10 7 m 2 /s 2 is the Earth's potential at the reference geoid while g denotes the Earth's 
acceleration (gravitational plus centrifugal), and where he is the clock's altitude above the reference geoid. 

Finally, we present the relation of the proper time read by the clock on the surface of the Earth at point C with 
respect to the TCB. Expressing drc/dt = (rfrc/^TCG) (dtTCG/dt) , with the help of Eq. © together with Eqs. ([!§]) 
and (|23|) . at the level of accuracy sufficient for GRAIL, we have: 



dr c _ 1 
dt ~ c 2 



i(v E + [n E xy c ]) 2 + C/ E (y c ) + N U b (r bE +y c ) + O(10" 17 ), (25) 



where the first term in the brackets, ve + [fi E x y c ] = ve + vco = vc, is the barycentric vel ocity of the DSN 
station. For details on the recommended relativistic formulation of GCRS consult Refs. [ll|, EH Il8| - Coordinate 
transformations (in particular, transformations involving topocentric coordinates) are discussed extensively in the 
IERS Conventions 2 . 

D. Relativistic timekeeping in the Solar System 

Spacecraft radio science observations are clock and frequency measurements made at Earth stations 18]. For this 
purpose, the time coordinate called Terrestrial Time (TT) is defined. TT is related to TCG linearly by definition: 

MTT =1-L G , (26) 



where Lq = 6.969290134 x 10~ 10 by definition. This definition accounts for the secular term due to the Earth's 
potential when converting between TCG and the time measured by an idealized clock on the Earth geoid [Tl] - [l3i [l8| . 
Using Eq. p3|) . we also have 



dr c drc dt TC G , , r 1 



dtxT dtTCG dtTT G C 2 



2 



Vco + CMy c ) +0(0-*). (27) 



On the other hand, equations of motion in the Solar System arc often evaluated using another defined time scale, 
TDB. TDB time (£tdb) is also related to TCB time t linearly: 

1.. w> 

where Lb = 1.550519768 x 10 -8 by definition, accounting for all secular terms due to the solar gravitational field, the 
Earth's orbital velocity, and the Earth potential on the geoid. 

The relationship between TCB and TCG is nonlinear; these are the coordinate times of two coordinate systems 
related to one another by the space-time transformations ([8]) -([9]) and their inverses (p~0|) (TTTT) . 

The relationship between TT and TDB, therefore, is also nonlinear. The difference is dominated by an annual 
periodic term with an amplitude of ~ 1.6 x 10 -3 s. The definition of TT and TDB ensures the absence of a significant 
linear term. 

For accurate computations in the SSB reference frame, observed times of transmission and reception need to be 
converted from TT to TDB. 

E. Coordinate reference frames in the vicinity of the Moon 

In the vicinity of the Moon, once again we consider two coordinate systems. The lunicentric LCRS is a coordinate 
system used, for instance, to represent lunar orbits. To describe experiments carried out on board the GRAIL 
spacecraft, we use the SCRS. 



2 All software, technical specification and other relevant materials associated with the IERS Conventions (2010) can be found at 
http : //www . ier s . org/ 



9 



1. Lunar Coordinate Reference System (LCRS) 



In complete analogy to the formulation of the GCRS (discussed in Sec. Ill C 1|) and similarly to the approach 
advocated in Ref. Q, the metric tensor of the LCRS may be obtained by transforming the metric ([3]) and the 
potentials (0])-© of the BCRS using the coordinate transformations given by Eqs. where it is sufficient to 

consider only the monopole contribution to U£ xt from all the bodies of the Solar System excluding the Moon: 

= £ ^ + 0(c-% a M = - VC& = -£ GM b ^ + 0{c-% (29) 



b^M 0IV1 b=£a 



r bM 



where summation is performed over all the bodies excluding the Moon (b =/= M) and the contributions due to the 
higher multipole moments of the mass distributions within the bodies are neglected. 

Applying the coordinate transformations given by Eqs. ©-© together with the external potential and acceleration 
(|29[) . one can derive the metric tensor g^ n of the non- rotating lunar coordinate reference system (LCRS). Denoting 
the coordinates of the LCRS as {y^} = (^M'^m)' this tensor may be presented in the following form (which is 
identical to the expressions in Sec. Ill CT1 after making the substitution E —> M): 

9m = l-^w M + ^w 2 M + 0(c- 6 ) 7 C = -7aA^M + 0(c- 5 ), g™ =i a p+i a ^w M + O(c- 4 ) 7 (30) 

where the scalar and vector potentials wm and w M are given as: 

w M = U M +ut\ dal + 0(c- 4 ), (31) 



G 



WM—rylyMxSMl+On, (32) 



where Sm in Eq. (|32|) is the Moon's angular momentum. Similarly to the GCRS, the scalar potential u>m (|3Tj) is a 
linear superposition of the proper gravitational potential of the Moon (with Rom being the Moon's radius, Mm its 
mass, while and Sfl are the Moon's spherical harmonic coefficients): 

oc +. 



u- 



M 



+ yy( ^Yp ek ( CO s evc^ cos w + sf k sin &</>)} + o( c - 4 ) (33) 



plus tidal contributions u^ dal produced by all the Solar System bodies (excluding the Moon itself, b ^ M) evaluated 
at the origin of the LCRF. For the GRAIL 's accuracy, it is sufficient to keep only the Newtonian contribution to the 
tidal potential produced by the external bodies which can be presented as 

< dal - E ( U ^ r bM + y M ) - U b (r bM ) - y M • VU b (r bM )) ]T ^ (s(n &M • y M ) 2 - y M ) + <D(y M , c" 2 ), (34) 

where n^M is a unit spatial vector in the body-Moon direction, i.e., n^M = ^bM/\TbM.\, where r^M is the vector 
connecting body b with the Moon. Note that the relativistic tidal contributions of 1 /c 2 order that are due to external 
potentials have a magnitude of 10 -16 when compared to Um and, thus, they were omitted in Eq. (|34[) . In addition, 
we present only the largest term in the tidal potential of the order of ~ y^; however, using the explicit form of the 
tidal potential Eq. (j3"4"|) . one easily evaluate this expression to any order needed for a particular problem. 

At the same time, we must account for deformations of the elastic Moon, expressed in the form of corrections AC^J 
and AS^, to the lunar spherical coefficients, due to the tidal potential of body b, located at lunicentric spherical 
coordinates (r bM , 4>bM, ObM.) [H,[ll : 



AC ik \_ M M b (RomY +L Ki + 2W-kW Jcob*W. 
AS ek j " 4A * Mm \T^J VTm)T" tk{ bM) IsinA&M/- (35) 



The lunar Love number k% — 0.025 [26l | thus introduces a significant time-dependent contribution to the spherical 
harmonic coefficients Cj^. and S^, which must be written as the sums 

C% = C%° + AC% and S™ = S™ + AS%, (36) 



where we used Cji° and SfP to denote the constant part of the lunar spherical harmonic coefficients. 
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2. Lunar Coordinate Time (TCL) 

There are several different time coordinates to be considered for GRAIL. In addition to the terrestrial time scales 
defined in Sec. Ill Dl GRAIL also relies on the timing events reported at proper times measured by clocks on board 
the lunar orbiters. Thus, one would need to introduce a realization of lunar coordinate time (TCL) and a spacecraft 
proper time (ST). 

The lunicentric orbits of the GRAIL spacecraft are coupled to the orbit of the Moon mostly through the difference 
between the acceleration of the probe and that of the Moon due to the gravitational pull of the Earth and the Sun 
(the Earth's and the Sun's tidal terms). This coupling is weak because the Earth and Sun tides are, respectively, just 
2.5 x 10 -7 and 4.7 x 10~ 8 times the monopole acceleration due to the Moon. Relativistic perturbations containing the 
mass of the Moon are small (~ 4.6 x 10~ n m/s 2 ) to the point that they are not measurable, being easily absorbed 
into the much larger non-gravitational perturbations (for instance, solar radiation pressure) . Should we conclude that 
general relativity does not matter in the computation of the the lunicentric orbit of the spacecraft? The answer is 
negative, but the main relativistic effect does not appear in the equation of motion. 

According to Eq. ©, the differential equation that gives the local proper time £tcl at the origin of the LCRS as 
it relates to the barycentric time TCB t is 



dt 



TCL 



dt 



= 1 - 



1 



V M 

2 



b^M 



bMJ 



) +0{cr% 



(37) 



where the Moon's barycentric velocity vm and position zm can be computed from the EIH equations ©. Equation ([57)1 
establishes the relationship between the TCL (£tcl) an d TCB {€) time scales. Truncated to the first post-Newtonian 
(1PN) order (we put the clock at the origin of its proper reference system, y M = and drop on the right hand side 
0(c~ 4 ) terms that are in principle known, but certainly not needed for our purposes), it is given by a differential 
equation 



^TCL _ 1 _ J_ 

dt ~ c 2 



2 



J2 GM b 



b^M 



TbM 



+ 0(c~ 4 ) « 1 - 1.48 x 10" 



(38) 



which can be solved by a quadrature formula once the orbits of the Moon, the Sun and the other planets are known. 



3. Satellite Coordinate Reference System (SCRS) 

To determine the metric tensor for the satellite coordinate reference system, the SCRS, we perform the coordinate 
transformation given by Eqs. where the "external" potential and acceleration determined by the potential 

w M given by Eq. (|3"Tj) , taken at the lunicentric position y A of the spacecraft GRAIL- A are (the equations are identical 
for GRAIL-B, except for the substitution A — > B): 

UL = wm(Ya) = ^m(Ya) + J2 ( U b( r bM + Ya) - U b (r bM ) - y A • V(7 6 (r 6M )) + 0(c' 2 ), (39) 

where y A is the solution of the equations of motion of the GRAIL-A spacecraft in the lunicentric frame. This equation 
can be obtained from equations of geodesies and the metric tensor of the LCRS (|30p with relativistic gravitational 
potentials given by (j3Tj) - (j3"2")l . (|3"3"|). and (j3"9")l . Including all the terms of the order of ~ 10~ 12 m/s 2 and larger, the 
equation of spacecraft motion of the GRAIL spacecraft in the LCRS takes the form: 

b^M 

GM M i 4(7 M- 

„2„2 



a A0 = -Vt/ M (y A ) - ]T (vU b (r bU + y A ) - VU b (i 

/4GMm ^ _ y2 + 4(nA _ VA())va \ + &NG + O ( 10 -13 m / g 2) (40) 

cy A v y A > 



where a A o = <i 2 y A /<ii 2 rcL and v A o = dy A /dtTCh are the lunicentric acceleration and velocity of the spacecraft in a 
non-rotating LCRS, also n A = y A /y A is the unit vector in the direction of the spacecraft, y A = |y A | and a^c is the 
contribution of non-gravitational forces affecting the motion of a spacecraft (e.g., solar radiation pressure, thermal 
imbalance, outgassing, etc). 

The first term in Eq. (|4T)]) is the contribution of the lunar gravity potential t/M(yA)i given by Eq. (|3"3")l . Note 
that the potential L/mCya) mus t be treated as the potential of an extended body and include a multipolar expansion 
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with sufficient accuracy. The second term in Eq. (|40)) is due to the tidal potential at the location of the spacecraft 
produced by external bodies wjj dal (mostly the Earth and the Sun) which is given by Eq. (f3~4"|). To reach GRAIL 's 
accuracy requirement of ~ 10~ 12 m/s 2 , one would have to account for several terms in the expansion beyond the 
second order one ~ given in Eq. In fact, terms up to ~ y^ in the tidal potential u^ dal are needed. The 

group of terms on the second line of Eq. ([40]) is the relativistic Schwarzschild perturbation due to the spherically 
symmetrical component of the Moon's gravitational field. The first two terms in this group are of the order of 
~ 1.86 x 1CP 10 m/s 2 and ~ 4.63 x 1CP 11 m/s 2 , respectively. These are large enough to be in the equations of motion. 
Given the nearly-circular orbit of the GRAIL spacecraft, the magnitude of the last term in this group is reduced by 
the orbital eccentricity, which is e A ~ 0.018. This fact reduces the contribution of this term by nearly two orders 
of magnitude when compared to the first two terms, making it barely observable with GRAIL. Note that the lunar 
angular momentum Sm present in the relativistic vector gravity potential Eq. (|32j) produces contribution to Eq. (|40j) 
of the order of ~ 10~ 14 m/s 2 , which makes it negligible for GRAIL. 

A a result, the metric tensor g^ an representing space-time in the coordinates {y A } = ($a> ^a) °f the proper non- 
rotating spacecraft coordinate reference system (SCRS) may be given in the following form: 



9oo 



= l-- WA + 0{c- b ), 
& 



9o a 



0(c-% 



9% 



7a/3+7a/3^^ A + 0(c 4 ), 



(41) 



where w A is the tidal contribution produced by the Moon on the world-line of the spacecraft: 

GM M /, 



w A 



2yi 



|3(n A .y A ) 2 -ylj+0(yi,c- 4 ). 



(42) 



We can also determine the differential equation that relates the rate of the spacecraft proper t a time, as measured 
by an on-board clock in lunar orbit, to the time in LCRS, £tcl: 



1 



2 

+ U u (y A ) + E ( U b( r bU + y A ) - U b (r bM ) - y A ■ VC/ h (r bM )) + (a M ■ y A )] + 0( C - 4 ), (43) 



b^M 



where a M is the barycentric acceleration of the Moon, Eq. (|29|) . see Ref. [151 ]. 

As a result, we can establish the rate of the spacecraft proper time with respect to the time of the BCRS, t = £tdb: 



dr A 
dt 



1 



1 



o(VM + v A0 

r bM + ) +0(c~ 4 ). 

b^M 



(44) 



This result summarizes the relationship of the proper time of an on-board clock in lunar orbit r A and the TDB. 



4- Transformation of gravitational potentials 

To complete the description of the LCRS, we present the transformation rules for the relativistic gravitational 
potentials. In Ref. we obtained the structure of the metric tensors corresponding to the local space-times in the 
reference frames relevant for GRAIL, expressed in terms of harmonic gauge potentials w and w. We also derived 
the rules for transforming relativistic gravitational potentials, the coordinate transformations between the frames 
and resulting relativistic equations of motion. Applying these results to GRAIL, we see that the scalar and vector 
gravitational potentials of the Moon wm ( v m) an d w M(yiv[) ( as measured at the LCRS) relate to those measured in 
the coordinates of the BCRS wm^bcrs) an d w M( r BCRs) as 

/ 2i> 2 \ 4 
wm(Ym) = (! + -^JwM(r B cRs) + ^2 ( v m • w M (r BCRS )) +0(c~ 4 ), (45) 

W M (Ym ) = W M ( r BCRS )- y M»M ( r BCRS ) + ° ) ■ ( 46 ) 

We estimate the magnitude of w M given by Eq. (|32|) as ~ 2.4 x 10 6 m 3 /s 3 . This results in a value of ~ 3.2 x 
10 -6 m 2 /s 2 for the third term in (1431) . which is four orders of magnitude too small compared to the second term in 
that expression, the scalar potential of the Moon multiplied by 2v^ A /c 2 that was evaluated to be ~ 5.5 x 10~ 2 m 2 /s 2 . 
Thus, to determine the relationship between the LCRS-defined mass of the Moon and its barycentrically defined mass, 
we must multiply the latter by the factor (1 + 2vf /l /c 2 ) « (1 + 2 x 10~ 8 ). Such a transformation results in a small, 
but observable effect. As far as the GRAIL's accuracy in concerned, contributions to other multipoles of the lunar 
gravity field are not sensitive to such a small correction factor. 



12 



As we know, GRAIL determines the lunar gravity field relying on the EIH equations of motion ([6]) for the bodies of 
the Solar System, including the Moon, the Earth, and the GRAIL spacecraft. For this, Eq. ^ must also includes the 
gravitational potential of the extended Moon (|33f and tidal potentials due to the Earth and the Sun (j34]) . Therefore, 
one would have to transform the resulting barycentrically defined gravitational potential of the Moon from coordinates 
of the BCRS to those of the proper lunicentric frame LCRS. A concern was that such a procedure may lead to some 
unwanted biases in the determination of the lunar gravity field. 

Our approach allows one to evaluate the general relativistic effects on the largest coefficients to the lunar gravity 
potential corresponding to this transformation. Substituting Eq. (|45[) into Eq. (|40l) , we essentially modify the equation 
of motion of the GRAIL spacecraft by accounting for the (1 + 2v 2 A /c 2 ) factor. Now, we can represent the barycentric 
velocity of the Moon as v M = v EM + v' M , where v EM is the barycentric velocity of the Earth- Moon barycenter and 
is the velocity of the Moon in the EMB frame. Therefore, v 2 ^ = v| M + 2v EM v' M cosD, where D = 9 — 9' 

is the the difference between the longitudes of the mean Moon and the mean Sun with a period of 29.531 days. The 
constant part in the barycentric velocity of the Moon vm may be easily absorbed in the determination of the lunar 
mass as a bias with magnitude of 2 x 10 -8 Mm- The variability in vm, if not properly removed in accordance with 
Eq. (|45"T) . may introduce an additional time-dependent bias with a magnitude of 2.2 x 10 _9 MmcosD, which can be 
removed in the data analysis. Finally, given the nearly circular obit of the GRAIL spacecraft around the Moon, the 
1 /c 2 terms in the second line of Eq. (|40|) can be seen and a modification of the Newtonian point-mass acceleration 
of the Moon a^J = GMm^a/v\- These terms are nearly constant, have combined mag nitude of - 1 x Kr 10 a^ and 
would be easily absorbed in the determination of the lunar mass as a small bias of 1 x 10~ 10 Mm- 

For spacecraft with lesser sensitivity, these corrections are irrelevant. However, at the micron-level sensitivity of 
the GRAIL mission, they become noticeable. It is, of course, possible to absorb small constant or periodic terms 
into constants such as Mm during data analysis, with no impact on mission objectives or the quality of the mission's 
results. Nonetheless, pursuing these small corrections is worthwhile, demonstrating that a spacecraft with GRAIL's 
sensitivity is already a practical instrument for relativistic geodesy in the lunar environment, and paving the way for 
future missions that will operate at even greater accuracy. 



F. Transformations of position vectors 

Equation (|9|) establishes the relationship between the coordinates of the local body-centric coordinate reference 
frame and the coordinates of the global BCRS (lfjj : 

r« = y a + c~ 2 { - |v Q (v a • y a ) - y a U^ t + [u a x y a ] + ±a Q y 2 - y a (y a ■ aj} + (D(c~ 4 ). (47) 

Considering the anticipated accuracy of the GRAIL experiment, we can simplify Eq. (l47l) by noting that the last 
three terms in this expression are much smaller than needed for GRAIL. Indeed, we can evaluate the magnitude of 
the third term 15] as [u> a x yj ~ y a G-M©i>EA£/AU 2 , where At is the signal propagation time, Mq is the mass of the 
Sun and AU ~ 1.5 x 10 11 m is the astronomical unit. The Moon-Earth radio-signal propagation time is At ~ 1.3 s. 
However, even for At = 10 3 s, we have [u a x y ] < 2 x 10~ A y a U£ xt ] therefore, the third term within the curly braces 
in Eq. (|47p is negligible for GRAIL. The last two terms within the curly braces in Eq. (|4"?| are dependent on the 
acceleration of the planet center, and are ignored for the reasons that we discuss at the end of this subsection. 

The required space-time transformations relate the position of the ground-based antenna and that of the spacecraft. 
Using superscript indices to indicate explicitly the dependence on the various time scales TT, TDB, etc., the terrestrial 
(geocentric) coordinates yJ T of the antenna must be transformed into TDB-compatible (barycentric) coordinates 
r c DB = z c — z e + 0(c ~ 4 ), where zq being the barycentric position of the DSN station zc = ze + Yq. Similarly to 
the approach developed in Ref. [27] for the BepiColombo mission, this transformation is expressed by 

r ™B = y TT^ _ 1 £ Ub{rm ^ _ 1 (v TDB . y TT) ^ (4g) 
C 6^E ° 

where X^e ^( x e) is the Newtonian gravitational potential due to bodies other than the Earth at the geocenter, xe 
is the SSB position and v™ B = dz-^/dt is the SSB velocity of the Earth (as defined in the paragraph before Eq. ((6|)). 

The time coordinate must also be changed consistently together with the spatial coordinates. The effect of this 
change on velocities is given by: 



TDB ,,TDB 



TT 
'CO 



1 - i £ U b (r m) ) - -L( V TDB . v TT )v T DB ««0 (49) 



b^E 



2c 2 



dt 
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where vJJ = c?yJ T /cLtq. Note that Eq. (|49|) contains the factor drc/dt, identical to Eq. (fl9l) . that deals with time 
transformation: tq is the local time for a ground-based antenna, that is, TT, and t is the corresponding TDB time. 
Similar to Eq. (|48|) . the lunicentric coordinates of the orbiter yA CL are transformed into BCRS coordinates r™ B = 
C(c -4 ), where za being the barycentric position of the spacecraft za = zr 



za ~ z M 



„TDB 



yl CL ' 



(l-^E^M))-^(v^ DB -yr 



TCLA TDB 



b^M 



J 



M 



(50) 



where X^b^M ^( x m) is the gravitational potential due to bodies other than the Moon at the Moon's barycenter, xm 
is the Moon's position in SSB coordinates and v™ B = dzu/dt is the Moon's SSB velocity. 
The corresponding velocity transformation is given by 



TDB 



TDB 
'M 



r TCL 
'AO 



TDB 
M 



TCL\ TDB 
A A0 J 



M 



b^M 



dt 



TCL 



dt 



(51) 



The relations for the other spacecraft is obtained by 



with v^cl = dy1 CL /dt T CL and c^tclMtdb given by Eq. {38] 
replacing A — > B. 

Note that in all the coordinate transformations presented in this section, we neglected terms that contain the SSB 
acceleration of the planet center, as these have an additional small parameter (R^/zb), where Rb is distance from the 
planetary barycenter and Zb is the distance of the planet from the SSB. Even for the Earth- Moon distance, these 
acceleration-dependent terms are at most of the order of 10~ 3 compared to the other 1/c 2 terms, both velocity- and 
Newtonian potential-dependent ones, making the acceleration-dependent terms negligible for the results above. 



III. FORMING Ka-BAND RANGE (KBR) OBSERVABLES FOR GRAIL 



One can demonstrate that in the post-Minkowskian approximation, appropriate for most Solar System experiments 
including GRAIL, as seen from BCRS, the phase of an electromagnetic wave that is passing by a gravitating body 
with the mass M can be presented as (see a detailed derivation in Appendix IA 21 with the result given by Eq. (|A22|) ): 



ip(t, x) = ko ( ct - Ra 



2GM 



■In 



Ri 



TA 



Ra 



0(G 2 ,c~ 



(52) 



where k m = k (l,k) is a constant null vector directed along the trajectory of propagation of the unperturbed 
electromagnetic wave such that ^ mn k m k n — 0, also k° = lo/c where tu is the constant frequency of the unperturbed 
wave. We also use the following notations: 



k = — , R A = x x A , 
Ka 



Ra = |Ra 



also 



r = x — z, 



r=r r A =x A -z, r A = r A 



(53) 



where z is the time-dependent spatial coordinate of the massive body. 

The phase <p of an electromagnetic wave that was emitted at the point x™ = (ciAOi x Ao) an d received at the point 
x b = ( c ^b, x b) remains constant along the path of this wave [28l 123]. In particular, if A is an affine parameter along 
the wave's path, the derivative of the phase satisfies the equation 



dtp 
d\ 



dip dx™ 
dx m dX 



= K m K m = 0, 



(54) 



which suggests that cp[x m (X)] = const. In other words, along the signal's world-line the phase stays constant and equal 
to its initial value ip(t,x) = kodAo = wao^ao- 

Equating the values of the phase given by Eq. (|52")l at two points Aq and B as <£>(£ao, XAo) = v(*b>Xb), we can 
determine the gravitational delay of the signal moving through a particular space-time. Indeed, up to 0(c~ 3 ) the 
coordinate time transfer, which is defined as £b — ^ao — Tab, is given by [l5| : 



Tab — ts — tAo — 



Rai 



2GM 



In 



^ao + r B + Rab 



tao + r B - Rai 



+ 0{c-% 



(55) 



where the logarithmic term represents the Shapiro time delay. Also, y-q = xb — z and Rab = xb — xao- 

With these results, we can now formulate a relativistic model for the fundamental timing observables on GRAIL. 
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SSB 

FIG. 2: Timing events on GRAIL: Depicted (not to scale) are the world-lines of the GRAIL-A and GRAIL-B spacecraft with 
corresponding proper times ta and tb. Ka-band signals that are emitted from points Ao and -Bo are received at B and A, 
respectively. The times tAo,^BO,iA and is are the coordinate times at these points, as measured in the BCRS. 



A. The inter-spacecraft ranging observables 

Consider a clock with proper frequency /ao, located at moving point Ao, that emits a signal with frequency /ao 
at an instant of proper time tao measured on the world-line of the clock. This signal is received by the moving point 
B at an instant of the proper time tb taken at the world-line body B and the instantaneous phase of this signal is 
compared with the phase of the local oscillator with proper frequency /bo of the clock located at point B. 

The measurable quantity is the difference between the instantaneous phases of the two signals compared at point 
B. Instrumentally, at point B one measures the fractional difference dn ab in the number of cycles dn^ received from 
the clock at point Ao and the number of the locally generated cycles dn B . Mathematically, this quantity may be 
expressed at the point B at an instance of the proper time c£tb as: 

dn ab = dn B - dnf Q = /bo^tb - /aV t b, (56) 

where f® is the frequency of the oscillator A as detected at B. 

Assuming that the number of pulses sent from spacecraft A, dn ao, and received on spacecraft B, dn^ , are the 
same, or dn^ = driAo, we can express /a via its value at the proper time of emission on spacecraft A (see also 
Eq. ([U3]) below): 

/ao _ dn Ao dr AQ = dTAo , g7 . 
/ao dr B dnAo dr B ' 

Furthermore, the instantaneous difference of the number of cycles measured on spacecraft B, as given by Eq. (|56l) . 
takes the form: 

dn AB = /bo^b - /ao*~ao- (58) 

Eq. (|58|) is the difference in the number of cycles generated by the two oscillators during the given proper time intervals 
along the world-lines of the two clocks. 

We can now express Eq. ([55]) in terms of the coordinate time: 

dn A B = /bo d *B - /ao (^~) *A0- (59) 

Note that Eq. (|5"9")l cannot be integrated in general case if the two time variables £ao an d t B are treated as indepen- 
dent. However, in our case the points Aq and B are connected by a time-like geodesic, and therefore, the coordinate 
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times £ao and t B are connected by the light-time equation (|55l) that reads: 



1 



t B - t A0 = Tab(*ao,*b) = -|r B (*B) - r A (iAo)| + 



2GM 



In 



ta + r~B - Rab 



(60) 



We note that Eq. (|60p can be used to express either t A o as a function of t B or vice versa. Observables on GRAIL are 
time-stamped using the time of reception (that is, t B ). We therefore have more direct access to x A (ts) rather than 
xa (t\o), and the first term on the right hand side of Eq. (|60| gets modified by Sagnac correction terms (as observed 
in Ref. [30]) consistently to the order 1/c 3 : 



-Rab = c^ae 



(d A B • v A ) , d A 



^ (vi + (n AB • v A ) 2 - (d AB • a A )) + 0(c~ 3 ), 



(61) 



where dAB — x b(^b) — x a(^b) is the coordinate distance between A and B at the moment of reception at B (we 
have cZab = |dAB | and nAB = d AB /d AB ) 7 where va = va(^b) denotes the coordinate velocity of spacecraft A at that 
instant, and where as is the acceleration of A (in all the order 1/c 3 terms we can use quantities at t A Q or is). In this 
case, Eq. (|60| becomes 



Iab(*b) = — — 
c 



(dAB • Va) 



^ab / ? , \? / , n\ 2GM , 

— (v A + (n AB • v A ) 2 - (d AB • a A )) + — g- In 



^A + ^B + ^AE 



r A + rB 



1AB 



+ 0(c- 4 ), (62) 



where all quantities here are taken at the instant of reception t B . In the case of the GRAIL mission, when signal 
transmission between the two spacecraft is concerned the first term in Eq. (|61[) . which is of order 1/c, is ~6.67 fis. 
The second term in Eq. (|BT|) represents the Sagnac term of order 1/c 2 and can amount to ^3.67 ns at GRAIL's orbit 
around the Moon; the third Sagnac term, of order 1/c 3 , is ^0.02 ps (comparable to the lunar Shapiro term, which is 
^0.04 ps). Note that expression for Tba(^a) may be obtained from Eq. (|()2"j) by interchanging A «-> B, 



B. Dual One- Way Range (DOWR) observables on GRAIL 

To develop an analytical form for the DOWR observable, we note that Eq. (|60|) could be used to express either £ao as 
a function of t B or vice versa. As observables on GRAIL are time-stamped using the time of reception (that is, t B ), in 
the following we treat £ao as a function of t Bl i.e., t A o = ^ao(^b)- Furthermore, we can write Tab (^ ao > ) = T AB (t B ). 
This allow us to present Eq. ([59]) as 

dn AB = /bo(^-) b *b - / A0 (^) A ^( iB ~ ^ab^b))- (63) 
To integrate Eq. (|6"3"|) . we rely on Eq. (l4"4l and introduce a function mb(^b) that allows us to write dr B /dt as: 

^ = l-\u B (t B )+0(c- 4 ), where u B (t B ) = 4 + E — + °( c ^)- ( 64 ) 

dt C 2 ^ ThA 

b 

Using this definition of u B (t B ) given in Eq. (|64[) allows us to integrate the first term in Eq. (|63|) as 

f t » /b °(S)b* B = /b °(S)b^ B ~ *^ + ^bo(«b(*b)(*b - t° B ) - JJ u B (t' B )dt' B ) + 0{c-% (65) 

where t B is used to denote the (for now, arbitrary) start of the integration interval. 
Similarly, we have the following expression for the second term of Eq. (|63p : 

fAo(^) A d(t B -T AB (t B )) = fAo(^) A (t B -T AB (tB)-(t ( B -T AB (tB))) + 

+ \fAo(u A (t B )(t B -T AB (t B )~(t B -T AB (t B ))^j- r AB(B) UA (t B X) +0(c- 4 ). (66) 

Transforming from proper to coordinate frequencies, as f B = f B o(dT / dt) B and f A = f A o(dr/dt) A , and using all the 
results developed in this section, we can integrate Eq. (|63j) and present the result in terms of the phase difference as 

An AB (t B ) = ht B - f A (t B - T AB {t B )^j + e AB + Sn AB + ©(c" 4 ), (67) 
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where e AB = cab(*b>*b) is given by 

e A B(t%,t B ) = ^{/Bo(uB(i B )(<B-t B )-^ o Wb(4)< 



fAa(u A (t B )(t B -T AB (t B )~ {t B -T AB (t B ))') - / B AB(B) "A(i B )^ B )}+0(^ 4 ), (68) 



rt& — Tab(^b) 

and Stiab = <5?iab(£b) ^ s an integration constant determined by the initial conditions: 

Sn AB (t° B ) = -/b^ + /a(^-T ab (^)) +0(c~ 4 ), 



(69) 



The second observable Auba(^a) that deals with the signal propagation from Bo to A is derived in an analogous way. 

To formulate the relativistic model for the dual one-way range (DOWR) observables on GRAIL, we need the 
expressions derived above for the instantaneous phase differences measured at both spacecraft, jtab^b) and ^ba(^a); 
which are given by Eqs. (|67[) . together with the instantaneous delays measured at the points of signal reception at 
both spacecraft, T AB (t B ) and Tba(^a), as given by Eqs. (1BT1) and (|62)l . As a result, Eq. (|67f becomes: 

An AB (te) = (/B-/ A )i B +.fATAB(iB)+eAB+<5nAB+0(c- 4 ) 



(/e 



•/a)*b + /a 
/a 



2GM] 



/ 4b _^ 

V c c 
(d A B • v A ) , d AB 



M 



In 



ta + r B + d A 



r A + r B - d AB 

2 



+ YT ( V A + ( n AB • v A ) 2 - (d AB • a A )) ) + cab + 5n AB + 0{cr 4 ), (70) 



where Mm is the mass of the Moon. An expression for Atiba(*a) may be obtained from (I7D1 by interchanging A ■<-» B. 

The authors of Ref. [6] discuss the interpolation algorithm realized on GRAIL to synchronize the LGRS clocks on 
both spacecraft in coordinate time. Here we just assume that synchronization is achieved, so that t A = t B = t and 
t B = t . We now can form a quantity, that is called dual one-way range (DOWR): 



l A 



An AB {t) + An BA (t) 



/a + /b 

Substituting Eq. ([70]) . we have the following result for Rdowr: 

((/ava ~ /bVb) ■ d A B) d AB (Cfeag - Ja&a) ■ &kb) 



(71) 



^dowr (^) 



^AB 



c(/a + /b) 
dAB/A(vi + (n AB -v A ) 2 ) 



+ 



2c 2 



/a + /e 



/b(v| + (n A B ■ v B ) 



2GM\ 



M 



2c 2 
c(e A B 



£BA 



/a + /b 



/a + /b 

?(<Wb + ^BA) 
/A + /B 



In 



?"A + r B + d Ai 

rA + r B - d Al 



+ 0(c-% 



(72) 



where e AB and 5n AB depend on the choice of the start of the integration intervals, i.e., t A and t B . 

We now discuss each of the seven terms present in Eq. (|T2"|) and evaluate their magnitudes and relevance for GRAIL. 
To develop numerical estimates for the magnitude of the various terms that we consider, we use mission parameters 
that are provided in Table fl] 

The first term in Eq. ([72jl is the instantaneous Euclidean distance c^ab — 200 km (see Table between the two 
lunar or biters. 

The next three terms are the first (~ 1/c) and the second (1/c 2 ) order Sagnac effects. These terms are due to the 
fact that representing the observables only in terms of the received times t A and t B on the two spacecraft is equivalent 
to a rotation of the reference system. To evaluate the second term in Eq. ([721 . we use the identity 



(/ava - /bv b ) = -|(/a + /b)(v b - v A ) - \{f B - /a)(va + v B ). 
Given A/ AB = |/b - /a| ~ 10 3 Hz, we get: 

((/av a - /bv b ) • d AB ) _ (v AB • d AB ) //b - /a\ ((va + v B ) ■ d AB ) 



(73) 



c(/a + /b) 



2c 



fA + fr 



2c 
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TABLE I: Select parameters of the GRAIL mission (some taken from Ref. @j) and the Earth-Moon system, along with 
corresponding symbols and approximate formulae used in the text. 



Parameter 


Symbol(s) 


Values used 


GRAIL Mission 


Inter-spacecraft range 


j 

«AB 


200 km 


Inter-spacecraft range-rate 


I^AB = (llAB ' Vab) 


2 m/s 


Lunar altitude 




55 km 


Lunicentric velocity 


VAO = VAO ^ VBO 


1.65 km/s 


Relative spacecraft velocity 


VAB — VAd.AB/(RM + ho) 


185 m/s 


Lunicentric acceleration 


oao = |aAo| — |aBO 


1.53 m/s 2 


Relative spacecraft acceleration 


fflAB — aAdAB/ (-Rm + ha) 


0.17 m/s 2 


Ka-band frequency 


/a =s /b 


32 GHz 


Frequency difference 


A/ab = /b — /a 


~ 10 3 Hz 


Earth-Moon system 


Moon's geocentric velocity 




1 km/s 


EMB orbital velocity 




30 km/s 


DSN geocentric velocity 




465 m/s 


Earth mass parameter 


GMe 


3.98 x 10 14 m 3 /s 2 


Moon mass parameter 


GM m 


4.90 x 10 12 m 3 /s 2 


Earth radius 


Re 


6.371 x 10 6 m 


Moon radius 


Rm 


1.737 x 10 6 m 



= (-0.061423 + 3 x 1CT 7 ) m, (74) 

where vab =vb- v A . Therefore, the second term in Eq. (|74[) is less than 1 /im and it can be omitted. 

The third term in Eq. ([72)1 is the second order (~ 1/c 2 ) acceleration-dependent Sagnac effect. We evaluate this 
term in a manner similar to Eq. (1741) and obtain the magnitude: 

d AB ((/bSb - Aa A ) • d A B) , (a AB ■ d A B) , / /b - Ia\ (( a A + a B ) ■ d A B) 

" "AB T^> "AB 



2c 2 / A + /b 4c 2 V/a + / B ; 4c 2 

= (2 x 10~ 8 + 1.2 x 10~ 14 ) m. (75) 

Thus, the entire third term in Eq. ([72"]) may be safely omitted. 

The fourth term on the right-hand side of Eq. ([72"]) is the second order (~ 1/c 2 ) Sagnac effect. As a result this term 
may be evaluated as 

c^ab /a(v| + (n AB -v A ) 2 ) + /b(v| + (n AB ■ v B ) 2 ) d A B {. 2 , / s2 , 2 , / \2 



■ 



2c 2 / A + / B 4c 2 V 

d AB / /b - / 



v A + (n AB • v A ) 2 + v B + (n AB • v B ) 2 + 



+ 4c^ ( / A + / B ) (( VAB '( VB + VA )) + ( nAB - VAB )( nAB -( VB + VA ))) = (0-002 + 2 x 10- 13 ) m. (76) 



Thus, the first term in ([76]) must be kept in the model. One can further evaluate this term by representing the 
barycentirc velocities of the GRAIL twins as v A = vm + v A o and v B = vm + v B o, where vm is the barycentric 
velocity of the Moon and v A o and v B o are the lunicentric velocities of the two orbiters. By doing this, one can see 
that there will be three terms, each of which is important for the GRAIL model. The term ~ c? AB (i , m/c) 2 contributes 
up to 2 mm to the DOWR. The term ~ d AB (vMV A o/c 2 ) contributes up to 110 /xm to the DOWR, and the last term 
~ d AB (v A o/c) 2 , a ls° frequency-dependent, contributes up to 6 /im to this observable. Thus, each of these terms must 
be accounted for in the relativistic model of GRAIL observables. 

The fifth term in Eq. ([72]) is the Shapiro gravitational time delay. Assuming a spacecraft altitude /iq = 55 km, this 
term contributes (4GM]y[/c 2 )(c? AB /(r A + r B )) = 12 /jm to the DOWR and, thus, it may be accounted for in the range 
model in the following approximated form, keeping just the largest (12 /im) term: 



2GA/ M 



In 



r A + r B + d A i 



r A +r B - d Al 



4GM M d AB , 4GM M d% B 



r A + r B 3c 2 (r A + r B ) 3 



12 nm+ 1.3 x 10 _B m. (77) 



Concerning the sixth term in Eq. ([72]) . in Appendix ]Bl we show that, for the times-scales of signal propagation 
realized on GRAIL (d AB /c ~ 1 ms), this term is of the order of 1/c 4 and contributes less than 1 x 10~ 15 (t — to) 2 m/s 
to the DOWR. An acceleration error of this magnitude yields a range error of less than 1 /im over the course of 6 
hours, and it is thus completely negligible. 
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The last term in Eq. (|72l) is of 0(c -2 ). This term represents the phase ambiguity in the DOWR observable at to- 
A method dealing with this term was outlined in Ref. |6|. We denote this term as Sn = Sn (t ) and keep it in the 
model. 

As a result, Eq. (fT2"j) can be presented in the following simplified form: 

iW W = d AB {l - (VAB o nAB) + ^ (v| + (n AB • v A ) 2 + v| + (n AB • v B ) 2 ) + 2 f° MM A + 0(0.5 M m). (78) 
^ zc 4c z V / c (r A + 7*b) j 

Up to this point, we treated the start to of the integration interval in Eq. (joTp)) as arbitrary. We now see that 
after negligible contributions are omitted, the start of the integration interval enters Eq. (178[) only in the form of 
the definition of the phase ambiguity Sn(t ). As we indicated above, dealing with this term is discussed in Ref. [f|. 
Once the effects of this phase ambiguity are accounted for, our formulation of the instantaneous DOWR observable, 
in the form of Eq. (|T8"|) . becomes independent of the choice of the start of the integration interval, and thus to is truly 
arbitrary, even as we maintain an instantaneous range accuracy better than 1 /im, as needed for the GRAIL mission. 

From Fig. Q] we can see that the vectors Ra and Rb are given as 

R A = x E m + x M + y A and R B = x E m + x M + y B . (79) 
These vectors are measured simultaneously with the signal reception in TBD and are needed to compute Eq. (|78p. 



C. Dual One- Way Range-Rate (DOWRR) observables on GRAIL 

To develop an analytical form for the DOWRR observable, we use Eq. ([59]) to express it as 

• u \ dnAB f f dTB \ r [ dT A\dt A0 

Using the notation / B = /bo(^ t b / c^b) and /a = f A o(dr A / dt A ) for the coordinate frequencies of the two clocks, we 
can present Eq. (I5U|) as 

% B (t B )=/B-/A^. (81) 

Mb 

As with the DOWR, the second observable DOWRR deals with the signal propagation from Bo to A is derived in 
an analogous way. Similarly to Eq. (|8ip . at the time t A on the spacecraft A we have 

nBA(t A ) = /A-/B$ £ , (82) 

dt A 

where the time of signal's emission i B0 mav be presented as a function of signal reception t A as t B0 = ^bo(^a)- 

Following the procedure outlined in Ref. Q, we assume that the LGRS clocks on both spacecraft are synchronized, 
such that t A = tB = t. We now can form a quantity that is called dual one-way range-rate (DOWRR): 

n AB (t) + n BA (t) / f A {dt A0 /dt B ) + fB{dt B o/dt A )\ 

^dowrr(t)=C — — =C[1 — — . (83) 

/A + /B V JA + JB > 

The ratio of coordinate times dt A o/dtB (and similarly dtBo/dt A ) can be computed by differentiating the coordinate 
time transfer equation (I5fl|) for tB — t A o = Tab^ao^b) (and similarly for t A — £bo = Tba(^A; ^bo)) with respect to 
the reception time £b ■ This procedure was already performed in Appendix IA3[ resulting in Eq. (|A33|) for dt A0 /dts ■ 
From this equation, the ratio dtBo/dt A is obtained by interchanging Af>B. 

Substituting these results for dt A o/dtB and dtBo/dt A into Eq. ([53"f we obtain the following expression for Wdowrr : 

M , , 1 f ((/bv b - /ava) • v A b) , ((/Ba B - f A a A ) ■ d AB ) 1 , 

^dowrr(i) = (n A B ' V A b) < 1 7 ~? ~ + " 7 ~i ~ \ + 

c{ f A + IB JA + JB J 

1 f / A (n A B ' V A )(v A B • Va) + /b ( n AB ' V B )(v A B ' V B ) 



C 2 I /A + /B 

4G c f M {rA ^ B rB)2 ((°A • va) + (n B • v B )) + 0(c" 2 ). (84) 
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The first term in Eq. ([84)) is the first order (~ 1/c) Doppler term, which may be as high as 2 m/s (see Table HJ; it 
clearly must be kept in the model. 

The second term on the right hand side of Eq. ([84)) can be evaluated using Eq. ([73)) as 



((/ B v B - / A v A ) • v AB ) v| b // b -/a\ ((va + v b ) -v A b) _ ., , ..„_,„., 



c(/a + /b) 2c VA + / B y 2c 



(5.2 x 1(T & + 6 x 10~ iu ) m/s. (85) 



Therefore, the second term in Eq. ([55]) can be dropped, but the v AB /2c term must be kept in the model. 

The third term in Eq. (l84l) is the second order (~ 1/c 2 ) acceleration-dependent Sagnac effect. To evaluate this term, 
we note that the acceleration vectors of the spacecraft point in different directions due to the ^200 km separation 
between the two craft. The vector difference can be calculated as &ab — 0.17 m/s 2 . We evaluate this term in a 
manner similar to Eq. (|75l) to obtain a magnitude of 

(ten - / A a A ) • d AB ) _ (a AB • d AB ) , ( / B - f A \ ((a A + a B ) ■ d^ = (g x 1Q _ 5 + 4 x 1Q _ U) ^ (gg) 



c(/a + /b) 2c V/a + Zb/ 2c 

We see that the second term in Eq. (|86l) is less than the needed accuracy of 1 /im/s and it can be omitted; the 
(aAB • dAB)/2c term, however, must be kept in the model. 

The (1/c 2 ) term on the second line of (|84|) can be presented as 

1 f /a(iiab • v a )(v A b • v A ) + /B(n A B ' v B )(v AB • v B ) -I If 1 

" J — ; f = ( n AB ' V A Vab ' V A ) + (n AB • V B )(v AB • V B ) \ + 

C 2 I JA + Jb > 2c 2 L J 



1 (h-l. 



2c 2 \h + fj 



) {(n AB • v B )(v AB • v B ) - (n AB • v A )(v AB • v A )} = (2 x lO" 6 + 6 x 1(T 14 ) m/s. (87) 



Therefore, only the first of the two terms on the right-hand side of Eq. ([87)1 must be retained. 

As a result, given the strict formation configuration implemented on the GRAIL mission, the model for DOWRR 
on GRAIL given by Eq. ([84)) has the following form: 



UdowrrW = (n AB • vab) - ^{ v ab + ( a AB • d AB )j + ^2 j(nAB ' va)(vab ■ v A ) + (n A B • v B )(v AB ■ v B ) j 

AGM M d AB 



(n A • v A ) + (n B • v B ) + °( - 1 (88) 



Equation ([88)) represents the instantaneous DOWRR observable for the GRAIL mission, developed to a level of 
accuracy better than 1 //m/s. One can verify that the result given in Eq. ([55)) may be obtained directly from Eq. ([75)) 
by simply differentiating Eq. (|78D with respect to time and retaining terms to the appropriate order. 



IV. CONCLUSIONS AND RECOMMENDATIONS 



We considered the formulation of a relativistic model for the observables of the GRAIL mission. We addressed some 
practical aspects of implementing the relevant computations. We derived an analytic expression that characterizes 
the process of forming the Ka-band ranging observables of GRAIL and developed a model for the dual one-way 
range (DOWR) observable. We also briefly addressed the transformation of relativistic gravitational potentials. This 
material can be used to improve the accuracy of modeling of the GRAIL fundamental observables. 

We presented a hierarchy of relativistic coordinate reference frames that are needed to GRAIL. In this respect, 
we introduced the barycentric (BCRS), geocentric (GCRS), topocentric (TCRS), lunicentric (LCRS) and spacecraft 
(SCRS) coordinate reference systems, together with the structure of the corresponding metric tensors in each of these 
systems and the form of the proper relativistic gravitational potentials — all presented at the accuracy required for 
GRAIL. We advocate a definition for the LCRS with its proper time, which we call the TCL. We presented the rules 
for transforming time and position measurements between the reference frames involved. 

The formula given by Eq. ([75)) is the main result of this paper. It is derived for the first time at this high level 
of accuracy including the terms of the 1/c 2 order. The final expression ([75)) is relatively simple and easy to utilize 
in practice. The equations we provide for time and frequency transfers are accurate to the level of 1 /im when used 
to analyze GRAIL ranging data. Modeling the DOWR observable at this level of accuracy is the most important 
priority for the mission and must be taken into account for the science data analysis. 

Most of the relativistic computations for GRAIL are done implicitly and are based on the models and tools available 
within the framework of JPL's Multiple Interferometric Ranging Analysis and GPS Ensemble (MIRAGE) software 
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[H . General relativistic equations of motion form the "back-bone" of the entire suite of models in MIRAGE and rely 
on the formulation given in Ref. [la ]. To navigate the GRAIL spacecraft, the code transforms the proper time of each 
of the GRAIL spacecraft to the time based on the SSB frame and integrates the spacecraft's barycentric equations of 
motion. To determine the inter-spacecraft range, the code then iteratively solves the barycentric light-time equations 
in terms of instantaneous distance (by recomputing the transmitter's position bearing in mind the elapsed light-time) 
in the presence of the Shapiro term. The analytical closed- form solution for DOWR (|78|) is not only more elegant, it 
allows for direct investigation of the observables and possible error terms under various circumstances in data analysis. 

We also developed a similarly accurate formulation for the DOWRR observable. Equation (|88[) allows us to calculate 
the value of this observable with an accuracy that is significantly better than 1 /im/s. In the form presented, Eqs. (|78|) 
and ([55]) can be readily incorporated into computer code used to model orbits and radio-metric observables. All the 
quantities in these equations are directly computable once the numerical positions and velocities of the spacecraft and 
Solar System bodies are known. It is also relatively straightforward to compute partial derivatives of these equations 
with respect to the the GCRS, which facilitates their use in fast numerical integration codes and optimizing solvers. 
Lastly, although we presented ideal, noise- free solutions, one can add relevant noise sources, including those in Ref. @. 

In a practical sense, the small relativistic terms that we calculated are easily absorbed into constant and periodic 
ad-hoc biases that are introduced during data analysis, with no impact whatsoever on mission objectives or the quality 
of the mission's results. Yet the existence of these terms, and the fact that they are observable at the level of sensitivity 
of the GRAIL mission demonstrate that GRAIL is already a practical instrument for relativistic geodesy, especially 
after the mission is complete and all data sets for the primary and extended mission phases are assembled 4, 5]. For 
future spacecraft that operate at even greater accuracy accounting for these relativistic terms will be essential. 

Although this paper was aimed specifically at discussing the range and range-rate observables of the GRAIL mission, 
we note that the solutions presented here are also applicable to other, similar missions. Foremost comes to mind the 
GRACE with a mission design very similar to that of GRAIL. Indeed, the calculations presented here may help shed 
light on the origin of small residual terms that were seen in the GRACE range and range rate observables [31). We 
will further investigate this possibility with results to be reported elsewhere. Clearly, the model to be developed for 
the GRACE Follow-on mission [H, Hj] must include similar higher-order terms to reach the anticipated DOWR and 
DOWRR at the level of few nm and nm/s correspondingly. We will address these issues in a subsequent publications. 
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Appendix A: The phase of an electromagnetic signal in gravitational field 

In this Appendix, we present derivations of the formulas use for time and frequency transfer in the GRAIL experi- 
ment. The derivation presented in this appendix is based on material that can be found in standard textbooks such 
as Refs. [2{| and (24|. A general solution is presented to the problem of light propagation in a gravitational field in 
the linearized approximation. 

1. General-relativistic post-Minkowskian space-time 

To develop the solution to the equations of the general theory of relativity in the post-Minkowskian approximation, 
we introduce the post-Minkowskian decomposition of the metric tensor g mn as 

9mn = Jmn + K ln + C(G 2 ), (Al) 

where h mn denotes the post-Minkowskian perturbation of the Minkowski metric tensor 7 m „. Following (29j . we impose 
the harmonic gauge condition on the metric tensor g mn , given in the form 

d m (V=9g mn )=0, or d m h mn - \d n h = 0(G 2 ), (A2) 

where h — jkih kl + 0(G 2 ). In the first post-Minkowskian approximation of general relativity [20l . [29J, Einstein's 
equations R mn = lQirG / c 4 (T mn — ^g mn T^ take the following form in arbitrary harmonic coordinates {x m } = (ci,x): 

(J?^2 ~ y2 ) hmn = ^r- ( Tmn ~ h mnT ) + °i° 2 )^ ( A3 ) 

where T mn is the stress-energy tensor describing a body that deflects a light ray and T = ^mT + 0(G). In the 
linearized approximation and neglecting the higher multipole moments, this tensor is given as [281 ]: 

T mn (t,x.) = Mu m u"v/l-/3 2 (5 3 (x - z(i)) + 0(G), (A4) 
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where M is the rest mass of the body, z(t) its time-dependent spatial coordinate, f3 = c x dz/dt and <5 3 (x) is the 
three-dimensional Dirac delta function. The body's normalized four- velocity u m , such that u m u m = 1, is given by 

.0 „, a \ „.o ft va \ „.o_ dt _ 1 , ^,n\ „.«_ uQ / c 



u m = (ti°,u Q ) =u° 1,— , u" = — = - +Q(G), u a = ' = + Q(G). (A5) 

Note that in Eq. (|A4|) we have neglected the factor of y/—g. This is done because in the linearized approximation 
\J—g = 1 + 0(G) and the quadratic terms cx G 2 are irrelevant in T mn since the corresponding time-dependent terms 
of the second post-Minkowskian order are currently unobservable in measurements made in the Solar System. 
We can now write down the Green's function solution to Eqs. (|A3I) : 

4GM 



h mn (t, x) = J y/1 - /3' 2 (u /m u' n - i 7 m ") <5 3 (x' - z(t'))G(i, x; t', x')d 3 x'dt' + 0(G 2 ), (A6) 

where G(i,x;£',x') is the Green's function 

G(t,x;t' ) x / ) = G(t-t';x-x') = — , 1 n d(t-t' - -Ix-x'lY (A7) 

4tt |x — x I V c / 

Integrating Eq. (|A6j) . one obtains the post-Minkowskian metric tensor perturbation in terms of retarded Lienard- 
Wiechert tensor potentials [29l . [33 | : 

^(t,x) = 4 T" £ + 0(G), (A8) 

c z u m r m 

where r m = x m — z m (t rct ) — [c(t — i rct ),(x — z(t ret )]. In Eq. (IA8|) . all time-dependent quantities are taken at a 
retarded time f re t (defined by Eq. (|A10I) below): u m = u m (t TCt ) = c~ 1 dz m (t Ict )/dt TCt is the body's four- velocity, 
/3(i r ct) = c~ 1 dz(t le t) /dt re t is the body's coordinate velocity normalized to the speed of light c. 

In the solution to Eq. (|A3I) . given in terms of the retarded Lienard-Wiechert potentials, all quantities involved 
including the distance r m = x' m — z m (t lct ), the body's world-line z m (t Iot ) = [ct re t, z(i re t)], and the four-velocity 
u m (t TC t) are functions of the retarded time t TCt . It is known (see, for instance, [34j ) that the retarded time in the first 
post-Newtonian approximation may be found from the null-cone equation 

lrnn r m r n = lmn [x m - z m (t Ict )][x n - z n (t Ict )} = 0, (A9) 

suggesting that the retarded time t ret = t ret (t,yi) is established as a solution to the equation 

*r«t = *--|a;-*(t*rt)|. (A10) 
c 

Note that Eq. (|A10p has an analytic solution only in the case of uniform motion of the gravitating body along a 
straight line. 

2. Phase of the electromagnetic wave 

The phase of an electromagnetic wave is a scalar function that is invariant under a set of general coordinate 
transformations. In the geometric optics approximation, the phase is found as a solution to the eikonal equation 

g mn d mV d nV = 0, (All) 

with g mn — 7 m ™ — h mn + 0(G 2 ). Equation (|A11[) is a direct consequence of Maxwell's equations. Its solution 
describes the front of an electromagnetic wave propagating in curved space-time. The solution's geometric properties 
are defined by the metric tensor (IA1[) . where h mn (|A8I) is the solution of the linearized Einstein equations (|A3[) with 
stress-energy tensor (|A4I) . 

To solve Eq. (jAlip . we introduce a covector of the electromagnetic wavefront in curved space-time, K rn — d m tp. We 
use A to denote an affine parameter along the trajectory of a light ray being orthogonal to the wavefront ip. Vector 
K m = dx m /dX = g mn d n <fi is tangent to the light ray. Equation (|Allj) states that the vector K m simply is null or 
gmnK m K n = 0. Therefore, the light rays are null geodesies [28| described by 

^ = \d m9kl K k K l . (A12) 
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Since eikonal and light-ray equations, given by Eqs. (|A11[) and (|A12[) respectively, have equivalent physical content 
in the general theory of relativity, one can use either of them to study the properties of an electromagnetic wave. 
However, the eikonal equation offers a more straightforward way to study the propagation of a wave. To find a solution 
of Eq. (|A11[) . we expand the eikonal ip with respect to the gravitational constant G assuming that the unperturbed 
solution of Eq. (|A11I) is a plane wave. The expansion may be given as 

<p(t, x) = + / k m dx m + <p G (t, x) + 0(G 2 ), (A13) 

where tpg is an integration constant and k m = fe°(l,k) is a constant null vector (i.e., ~fmn.k m k n = 0) along the 
direction of propagation of the unperturbed electromagnetic wavefront. Furthermore, k° = lo/c where lo is the 
constant frequency of the unperturbed wave, and tpc is the perturbation of the eikonal to 0(G), which is yet to be 
determined. Substituting Eqs. (|A8[) and (|A12I) into (jAll \ and keeping only first order terms in G, we obtain an 
ordinary differential equation to determine ipQ\ 

fe.Wfe.^flS^+o^,, (A14) 
dX 2 cr u m r m 

which alternatively can be obtained as a first integral of the null geodesic equation (|A12p . We can now integrate 
(1A14[) while keeping in mind that dk m = and employing the exact relationship [34j : 

"' A ' d\n[k m r m }. (A15) 



ri I /ytTYL fsi qrtTFb 1st fit 

Urn' i^'m 1 rvmw 

Neglecting the body's acceleration (or du m — 0), a plane- wave solution of Eq. f|A14[) has the form 

2 C 1 Af 

p G (t,x) = -^-(k m u m )\n[k m r m ] , (A16) 

where all quantities on the right-hand side are taken at the retarded instant of time t re t in agreement with (|A10p . 
Therefore, we can now write the post-Minkowskian expansion for the phase of the electromagnetic wave as: 

<p(t,x) = <p +l k m dx m + ~^{k m u m ) In [k m r m ] + 0(G 2 ), (A17) 

which can be presented in the following form 

<p(t, x) = po + fc (ct - k • x + 2G 2 M1 7 (V ' k)/ | ln [k°(r - k • r)] ) + 0(G 2 ), (A18) 

where all the quantities in the last term are taken at the retarded time t rct defined by Eq. (lAlpp . 

Let us now consider signal propagation from a point (c£ai x a) to a point (ct, x). Then along the signal's path the 
phase (|A18[) will change according to 



V {t, x) = k U - R A + g^ 1 -/ v - 2 k )/ c ln ) + 0(G% (A19) 



c 2 y/l-v 2 /c 2 



k • r 



r A - k • r A 



where we used the following notations: k = Ra/-Ra, R-a = x X A, = |Ra| and r = x — z(i rct ), r = |r|. 
One can further simplify the argument of the logarithmic term in Eq. (|A19[) as 

r-k-r r A +r-R A 



ta - k • ta r A +r + R A 
Thus, Eq. (|A19I) takes the form: 



2GMl-(vk)/c. 



y/i ~V 2 /C 2 



r A + r + R a 1 



ip(t,x) = k [ct -R A 5 ■=. v = '' In — ^ \+0(G 2 ). (A21) 



r A + r - R A 



Effects of order of vG/c 3 are very small and one can neglect them in Eq. (|A21|) . As a result, the post-Minkowskian 
phase of the electromagnetic wave, with an accuracy appropriate for modern-day Solar System experiments [30, [3J] , 
can be presented as: 

r A +r + R A 



2GM . 

in 

-r A + r - R A . 

Along the signal's path the phase stays constant and equal to tp{t A , xa) = kod A . 



^(t,x) = Uct- R A ^iln [ rA t r + ^ A D +0(G 2 ,c- 3 ). (A22) 
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3. Coordinate gravitational time delay 



Consider the case of one-way signal transmission. Let A be the emitting station, with BCRS position x A (t), and B 
the receiving station, with position x B (t). Also, z(t) is the vector connecting the SSB to a gravitating body. We denote 
by f a the coordinate time at the instant of emission of a radio signal, and by t B the coordinate time at the instant of 
reception. We put r A = x A (i A ) - z(t A ), r B = x b (*b) - z(*b), also Rab = x B (* B ) - x a (*a), N A b = Rab/-Rab and 
i?AB = IRabI, r A = l r A I j r B = | r sl are the Euclidean norms of these vectors. 

We know that along the signal's path the phase stays constant. Thus, equating the eikonal of the wave given by 
Eq. (|A22|) at the two points A and B as (p(t A ,x A ) = </?(£b,x b ), we determine the gravitational delay of the signal 
mo ving through a stationary space-time. Indeed, up to 0(c -3 ), the coordinate time transfer Tab = £b — ^A is given 
by EE}: 



R AB 2GM 
Tab =t B -t A = 1 — In 



rA + r B + Ray 



ta + tb- Rai 



(A23) 



where the logarithmic term represents the Shapiro time delay. 

The ratio of coordinate times dtAo/dt B (and similarly dt B o/dt A ) can be determined directly by differentiating the 
coordinate time transfer equation (|60p for t B — £ao = T AB (t A Q,t B ) (and similarly for t A — t B o = T BA (t A ,t B o)) with 
respect to the reception time t B . In other words, we must evaluate 

d , , d rl. , , , . , 2GM , / ta + r B + Rab \} , . „ ,% 

-r-{t B - t A0 = -T- - re iB - r A (t A o + — 3- In , „ ) • (A24) 

dt B dt B lc c a V ta + r B — Rab ' > 

While performing the differentiation, we must account for the fact that the coordinate distance between Ao and B 
depends on both the times of emission and reception, i.e. Rab — |i"b(£b) — i"a(£ao)|- For instance, we have 

dRAB „ / dt A0 \ 

Nab- vb-va— I. (A25) 



dt B V dt 

For notational convenience, we henceforth denote the vectors by ta = xa(^a) — z(2a) and r B = x B (i B ) — z(t B ). 
As before, we have r A — \r A \ and r B = |r B |, as well as the coordinate velocities va = xa(^a) — z(£a) and vb = 
x B (i B ) — z(t B ). Performing the differentiation in Eq. (|A24j) to order 1/c 3 , the last factor in Eq. (IClj) dt A o/dt B 
can be given by the ratio (in agreement with the recent work on the ACES [30( and RadioAstron [35[ missions; for 
convenience, we adopt notations similar to those introduced in Ref. |30|): 



dtAO _ 9b 
dt B q Aa ' 



(A26) 



where q A o and q B are derived from 



i 1 /tvj x 4GM (r A + r B )(N AB • v A ) + Rab{t A ■ v A )/r A , > 07 , 

ffA o = 1 - -(Nab ■ v A ) - — ^ + ^ _ ^- , (A27) 

1 4GM (ta + r B )(N AB ■ vb) - Rab (£b ■ v B )/r B 
«b = 1 - -(Nab ■ v B ) - . (A28) 

In an experiment, the position of the transmitter A may be recorded at the time of reception t B rather than 
at the time of emission £aoj i- e - we may have more direct access to yL A {t B ) rather than xa(£ao)> and the formulae 
()Al7)) - (|Al8t get modified by Sa gnac correction terms consistently to the order 1/c : 



<7ao = 1 - — ( n AB • va) - \ fv| - (iab • v A ) 2 - (a A • Af 
c c z V 

+ ^j{( v A - ( n AB • v A ) 2 )(n A B • va) + d AB {i{eL A ■ va) - (n AB ■ a A )(n A B • v A ) - (a A • d AB )) } 
4GM (r A -I- r B )(n A B • v A ) + c^ab (i"a ■ v A )/r A _ 4 



(r A + r B ) 2 - 4 



+ C(c- 4 ), (A29) 



AB 



«b = 1 - -(nAB ■ v B ) - -4 ( (v A • v B ) - (n A B • v A )(n A B ■ v B ' 



2^3 {( V A ~ ( n AB • v A ) 2 )(n A B • v B ) + <^ab ^( a A • v B ) - (n A B ■ a A )(n A B • v B )) | 
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AGM (r A + r B )(n AB • v B ) - d AB (r B ■ v B )/r B + ( c -4) ( A30 ) 



(r A + r B ) 2 -d 2 



AB 



For the GRAIL spacecraft in lunar orbit, in the lunicentric reference frame the first, 1/c term in Eqs. (|A29[) - ([A30[) 
is ~ 10~ 4 in the bary centric frame and ~ 5.5 x 10~ 6 in the lunicentric frame. The second term is the 1/c 2 Sagnac 
term, which is ~ 10" 8 for the BCRS and ~ 10 -11 for the LCRS. The 1/c 3 terms are ~8x 10 -17 for the Sagnac term 
and ~ 1.4 x 1CP 16 for the Shapiro term. Note that the result given in Eqs. (|A29[) - (|A30[) has been obtained assuming 
that the field of the Earth is spherically symmetric. Indeed, the J2-terms in the factor (<7a/<Zb) do not exceed 10~ 15 . 

As a result, the expression dt A0 /dt B from Eq. (|A26|) has the form 

dtAQ qB , 1/ N 1 (, \ i / a \\ 

-7— = = 1 (riAB • vab) j ( v a • v AB ) + (a A ■ d A Bj - 

dt B qAo c c A V / 

- { 2 ( V A • vab)Kb • va) + (v A - (n A B • v A ) 2 )(n A B • v A b) + 

+ d AB (2(a A • va) + 2 (nAB ■ a A )(n A B ■ v A ) - (a A • v A b) - (nAB ■ a A )(n A B • v A b) - (aA • d AB )J } - 
AGM (r A + r B )(n AB • v A b) - d AB ((r B • v B )/r B ) + (r A ■ v A )/r A )) 



(r A + r B ) 2 -d 2 



+ 0(c- s ). (A31) 



AB 



We can evaluate the magnitude of each of the seven terms in Eq. (IA31[) using the values from Table |U We will use 
these basic formation parameters to evaluate the terms in (j8~4")l . keeping only those that contribute to the range-rate 
model more than Aw^owrr — 1 A*m/s. The error terms in these expressions must be less than Av^ owlT /c = 3 x 10 -15 . 

Keeping these numbers in mind we see that the (1/c) term in (|A31[) is of the order 7 x 10 -9 and must be kept in the 
model. The two 1/c 2 terms were evaluated to be of the order of ((va • vab) + ( a A • d AB ))/c 2 ~ 6 x 10 -11 + 3 x 10 -12 . 
Thus, both of these terms must be in the model. 

Among the 1/c 3 terms, the first one is ~ 6 x 10~ 15 and must be kept. Next is the term that contains (nAB • vab), 
which was evaluated to be at most 3 x 10~ 17 and is too small to be in the model. The remaining 1/c 3 term, which is 
oc c?ab, is at most ~5x 10~ 16 and thus, this entire term may be neglected. 

Lastly, the Shapiro term in Eq. (|A31|) was evaluated to be 

AGM (n AB • v AB ) . AGM d AB ((n A ■ v A ) + (n B ■ v B )) _ 19 lfi 



(rA + r B ) c 3 {rA+r B ) 



8 x 10" 19 + 3 x 10" 15 . (A32) 



Therefore, we will keep only the second part of the Shapiro term in the model. 

As a result, the expression (|A31[) may be presented in the following simplified form: 



dtAO 



1 - -(nAB • v A b) - \ ( (v A ■ v AB ) + (a A • d AB 
c c z \ 



dt B 

- 1 ((n AB • v A )(v AB • v A )) + 4GM d AB ( {ua . Va) + (nB . Vb) \ + 0{5 x 10 -i 6)) (A33) 
c 6 V / c 3 (rA + r B y \ / 

where nA = ya/ta and n B = r B /7' B . This expression accounts for all the terms that may have magnitudes larger 
than 5 x 10~ 16 and are thus relevant to the GRAIL mission configuration. 

Appendix B: Evaluating the DOWR integral expression 

In Sec. IIIIBI we defined the quantities cab(^b) and eba^a), which were given in Eq. (1681) . We now construct the 
quantity 6ab(*b) + £ba(*a): 



cab 



{t B )+e B A{t A ) = -^{f bo (u B (t B )(t B -^)-u B {t A )(t A -T BA (t A )- (*a - r BA(ti))) 



*b— *a+T b a(*a) 

u B (t')dt' 

t° B -t° A +T BA (tl) 



+ /ao(«a(*a)(*a - *a) - «a(*b)(*b - T AB (t B ) - (t° B - T AB (t B ))) - 

B( * B) uaW)|+0(c- 4 ). (Bl) 

(tg) n 



*a— *b+Tab (*b) 
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To estimate the magnitude of this quantity, we can assume that the clocks A and B are perfectly synchronized. 
In reality, nothing is perfect and the GRAIL mission's design relies on a synchronization procedure (see Ref. @ for 
details) that, of course, leaves synchronization errors that percolate in the data analysis. One can study the impact 
of synchronization errors on the measurement accuracy, using the approach outlined here. However, presently we 
concerned with the evaluation of the magnitude of the error that may arise in the ideal case, if we were to drop the 
£ab and eba terms. Thus, we assume that t A = t B = t and t A = t B = t , so that the quantity given by Eq. (|Bip is 
reduced to 

Tab (^) 

e AB (t) + e BA (t) = \\f A0 (u A (t)(T AB (t)-T AB (t )) - [ u A (t')dt') + 
c L v JT AB (t«) ' 

+ f bo {u B (t){T BA (t)-T B A (t )) - / u B {t')dt') \+0{c- i ). (B2) 

V JT BA (t°) J } 

We series expand the integrals in Eq. (|B2I) around the start t° of the integration interal. For the integral term 
multiplied by /ao inside the second set of square brackets in Eq. (IB2|) . we have u A (t') = u A (t°)+u A (t°)(t' — t )+O(At 2 ) 
(where At = t — t°) , and obtain: 

fT AB (t) 

u A (t)(T AB (t) - T AB (t )) - / u A (t')dt' = 5 u A (to) (Tab (t ) - Tab (to)) + 0{At 3 ). (B3) 

•/Tab (to) 

The integral term in the first set of square brackets in Eq. (|B2[) . which is multiplied by /bo, can be evaluated in a 
similar way. Keeping just the leading terms (~ d AB /c) in Tab and Tba, we present, for instance, T AB (t) — Tab (to) — 
c _1 dAB(t - to) + 0{cr 2 ), so that Eq. jB2j} becomes: 

e A B(t)+e B A(t) - ^(f A ou A (t ) + f BO Mto))d 2 AB (t-t ) 2 + 0(c- 4 )=0(c- 4 ). (B4) 

Remembering the form of u B (t) in Eq. (|64|) . and denoting by /q the typical GRAIL radio frequency (/ao — /bo — /o, 
see Table[T|, we see that the term above is of the order of ~ (v A a A d\ B /2c A )(t — £o) 2 /o "C 7.7 x 10 -24 (t — to) 2 /o at 
most. This term is multiplied by ~ c/2/o as its contribution to the DOWR measurement is calculated in Eq. (|T2|) : 
the magnitude of this contribution is therefore less than 1 x 10~ 15 (t — to) 2 m/s , which is negligible for GRAIL. 

Appendix C: Relativistic frequency transfer between spacecraft and a DSN antenna 

The DOWR observable is defined in terms of the coordinate frequencies f A and f B of both signals generated on-board 
the two GRAIL spacecraft. However, these frequencies are measured by a ground-based DSN station. Specifically, in 
the case of the GRAIL spacecraft, while the spacecraft-to-ground transmission takes place using frequency bands that 
are different from the frequencies used for inter-spacecraft communication, the two frequencies are synthesized using 
the same timing source (USO) on board. Thus, measuring the frequency of the spacecraft-to-ground transmission is, 
in effect, a measurement of the inter-spacecraft frequency as well. 

Below, we simply assume that the spacecraft-to-ground transmission does serve as a means to measure f A and f B 
without going into further detail. We discuss how one can use these measurements, taken using the TT time coordinate, 
and transform them to the TDB time coordinate. We focus only on relativistic frequency transformations between 
the frames involved and will neglect frequency-dependent media effects, which are easy to reinstate when needed. 
Although communication between the spacecraft is conducted using Ka-band microwave signals and spacecraft-to- 
DSN is done relying on X-band signals, these frequencies are related by simple numeric factors. Therefore, in a slight 
abuse of notation, we use the same symbol for the inter-spacecraft and spacecraft-to-ground frequencies. 

1. One-way frequency transfer 

We consider the situation when a signal with frequency / is measured by an electronic counter whose register 
is incremented by 1 each time the magnitude of the signal changes from minus to plus. The number of cycles dn 
measured by this counter in the interval of proper time dr is then dn = fdr. Thus, the frequency transfer between 
transmitter (at point A) and receiver (at point C) requires the determination of the ratio /a//ao between the proper 
frequencies /ao transmitted by satellite (A) and f A on the ground (C). The infinitesimal proper time intervals dr A 
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and drc correspond to the infinitesimal number of cycles, dn, at the transmission and reception points A and C, so 
that diijl — dnA- Therefore, the one-way frequency shift during the transfer from A to C is 

/a dn^ drA drA (dr/dt) A dtA 



For (dr / dt)c / {dr / dt)A, we get 



/ao dr c dnA dr c (dT/dt) c dt c 
{dT/dt) A _l-c- 2 [U(r A ) 



(CI) 



(dr/dt) c 1 - c- 2 [U(r c ) + ±v 2 ] 



+ 0(c-% (C2) 



where U(rA) = J2b^b( r bM + Ya) an d ^( r c) = Sb^f>( r f>E + Yc) are the Newtonian potentials at the points A and 
C and vectors va and vc are the barycentric velocities of the spacecraft and the DSN station correspondingly. The 
factor given by Eq. (1C2[) consists of the Einstein gravitational red-shift and second-order Doppler effects, both of order 
1/c 2 . Therefore, Eq. (|C1|> becomes 

fg_ cWa _ 1 - c- 2 [U{va) + |vj] dtA . 4 . , > 

/ao rfr c 1 - c- 2 [[/(r c ) + iv 2 ] dt c + ^ j ' ^ 



The dtA I 'dtc term on the right-hand-side of Eq. (|C3|) is the ratio of coordinate periods of the same signal at A and 
C. It contains both the (~ 1/c) Doppler effect and the (~ 1/c 3 ) terms that we seek. To compute this factor at the 
accuracy of 3 x 10 -15 , it is sufficient to treat the Earth and Moon gravitational potentials as monopole potentials, 
i.e., the approximation U — GM/r and consequently, the use of Eq. (|44|) is sufficient. 

Similarly to the discussion in Sec lA3| the ratio dtA/ 'dtc can be computed by a direct differentiation of the coordinate 
time transfer Tac — tc ~ ^A with respect to the emission time tA- In Appendix IA 31 we computed all the relevant 
terms, accurate to the ~ 1/c 3 order. Using the time of reception tc (and expressing the time of emission as a function 
of the time of reception, tA = i A (ic))> the factor dtA/ dtc m Eq. (|C1[) can be given (similar to Eq. (|A26jl ) by 

dtA = qc (C4) 
dt c QA 

where qa and qc are then given to 0(c -3 ) from Eqs. (|A29|) - (|A30j) as 

9A = 1 - "Oac • v A ) - ^ ( v a - ( v a ■ n AC ) 2 - (a A • d A c)) + C(c" 3 ), (C5) 

qc = 1 - ^(n A c • v c ) - ^((v A ■ v c ) - (v A • n A c)(vc • n A c)) + C(c -3 ), (C6) 

where dAc — x c(^c) ~ x a(^c) is the coordinate distance between A and C at the moment of reception at C (we have 
dAc — |dAc| an< i n AC = dAc/^Ac)j where va(£c) denotes the coordinate velocity of the station A at that instant, 
and where ac is the acceleration of A. Therefore, the ratio qc/qA can be given as 

^ = 1 - Vac • vac) + ^ (Va - |v 2 c + |v| c - (a A • d AC )) + O^ 3 ), (C7) 

yA C C \ / 

where all the quantities involved are given at the time of reception tc- 
Therefore, for the one-way frequency transfer, we have 

ft l-cr*[U{vA) + \vl_ 



/ao 1 - c- 2 [U(r c ) + |v 2 ] g A ' 
with the ratio 9a/9c given to sufficient accuracy by Eq. (IC7[) . Up to C(c~ 3 ), this expression becomes: 



(C8) 



^ = l-i(n A c-VAc) + 4(5 v Ac+E[ C/b ( rhE +yc)-^(r 6 M + yA)] -(a A -d A c)) +0(5x 10- 14 ). (C9) 

J AO C C ^ / 

o 

Note that all the quantities in Eq. (|C9I) are taken at the reception time. This relation determines the frequency 
displacement due to the combined motion of the emitter and the receiver (intrinsic Doppler effect) and the difference 
in the gravitational potentials at the points of emission and reception of the signal (gravitational displacement of the 
frequency) . 

In the case of GRAIL, the one-way frequency transfer given by Eqs. (|C5|) - (|C6I) can be evaluated numerically as 
follows. The first-order Doppler effect is |(n A c • v A )/c| = 5.5 x 10~ 6 ; for the ground |(n A c • vc)/c| = 1.6 x 10 -6 . 
The second-order Doppler effect is v|/2c 2 = 3.4 x 1CU 10 for the satellite; v 2 -./2c 2 = 1.3 x 1CU 12 for the ground. The 
gravitational red-shift (Einstein) effect is given by U A /c 2 = U E (r A )/c 2 = 6.5 x 1CT 10 ; U c /c 2 = 6.9 x lCT 10 . The 
0(c -3 ) terms are less than 3.6 x 10 -14 for the spacecraft and 2.2 x 10~ 15 for the Earth station, thus, they are omitted. 
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2. Integrated Doppler Effect 



To measure the frequency of the received signal, DSN receivers count the number of cycles received from the 
spacecraft over intervals of time measured by high precision clocks located at the receiver station. 

We consider the integrated (one-way) Doppler effect that is at the basis of forming the Doppler observable for 
spacecraft such as GRAIL that are equipped with a precision on-board frequency source. Consider a clock (Fig. [2]) 
with proper frequency f A Q, located at point A (GRAIL-A), which emits light signals during the proper time interval 
(ta.1)TA2)- These signals are received with frequency f A by a ground-based DSN station, located at C (see Fig. [T]), 
over a proper time interval (tci,tc2)- The measured quantity over the interval (tci,tc2) is the number of cycles 
received at point C from the transmitter at point A. The number of cycles received at point C must equal the number 
of cycles emitted at point A. Hence, 



TC2 pTA2 

f A drc = / /ao^ta, 



(CIO) 



where the frequency /ao is the constant (with respect to ta) precision oscillator frequency of the transmitter, and 
may be moved outside the integral sign. The remaining integral can be transformed to the proper time re with the 
use of Eq. jCH]): 



dr A = 

with ratio drc /dt given by Eq. (|2~5|) as 

dr c 



' C2 dTA 



dr c 



1~C2 fC 



II dr c = I M-^-dt, 
n /ao Jt t /ao dt 



dt 



= 1 - — 



v 2 c + E u b( r bE + y c )l + o(io- 17 ). 



1,,2 

2 



(Cll) 



(C12) 



Also, in Eq. (|C9I) we expressed the ratio f A //ao a t the adopted accuracy using terms that are functions of t. Therefore, 



dr A 



d Ac(t2)-d A c(tl)) + 



hie - h v c~Yl u b(*bM+y A ) - (a A -d A c)l dt+0(c- 3 ), (C13) 



where we relied on the identity (riAc ■ vac) = c?ac- Using these intermediate results, we find that the total number 
of cycles N A , received during the count interval At = t% — ti, is given by: 



iVA 



T~C2 



f°drc = fAo(t 2 - h) - i/Ao(dAc(i 2 ) - d A c(h)) + 



+ -5T/ 



AO 



*2 r 



2" V AC ~ 



1 v c -^C/ b (r fc M+y A )-(aA-d A c) dt + 0(c^). (C14) 



Dividing the number of cycles N A by the count interval At yields the Doppler observable f A : 

fc _ £va 
/A At ' 



(C15) 



Therefore, the spacecraft frequency /a that is measured at a DSN receiver during the count interval of proper time 
At can be modeled as 



J A _ ^ _ 

/ao 



dAc{t 2 ) 



<M*i) +i 1 



cAt c 2 At j, 

where the range d A c is defined as d A c = l x c — X A 



2 V AC 



iv c - J2 u b(nu + y A ) - (a A • d AC )j dt + 0(c- 3 ), (C16) 



3. Transmitted spacecraft frequency, as seen from the BCRS 



In a 1-way Doppler transmission between a GRAIL spacecraft and a DSN station, the measured frequency received 
at the DSN station, f£, is related to the transmitted frequency /ao as given in Eq. (IC16I) . Our objective is to express 
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the frequency received at the DSN using the TDB time coordinate of the BCRS. To do this, we imagine a hypothetical 
receiver that is at rest with respect to the BCRS, and co-located with the DSN receiver at the moment of receiving 
the same signal. Using Eq. (|C9| we can relate the instantaneous frequency /a that is received by this hypothetical 
receiver to the transmitted frequency /ao'- 



h 



AO 



1 - ^(n A • v A ) + ^ (§v A + [U b (z b ) - U b (r bM + y A )] + (a A ■ d A )) + C(c~ 3 ), 



(C17) 



where J^b U b {z b ) is the gravitational potential at the origin of the SSB (see discussion after Eqs. (jU)-©) defined as 



5>£z ft = and ^ =G A/ b {l + ^(t|-^^)}+0( C - 4 ) 

b c^b 



(C18) 



also, z b is the barycentic position vector of the body b, v b = \i b \ is its barycentic velocity, and all times are measured 
in the TDB time. 

Analogously to the discussion in the Section IC 21 the Doppler observable /a that corresponds to the instantaneous 
frequency /a is established during the count interval of TDB time of At = t 2 — t\ and is described as 



h_ _ 1 d A (t 2 ) - d A (ti) 1 
/ao cAt c 2 At 



*2 r 



2 A 



Y [U b (z b ) - U b (r bM + y A )] + (a A • d A )l dt + 0{c^), (C19) 



where we used the fact that c?a = |^a| and c?a = (ha ■ va). 

Finally, with the help of Eqs. (jC16[) and (|C19[) . we obtain the following relation between two Doppler observables- 
the spacecraft's frequency /a reported at the BCRS and the same frequency / A as measured at the DSN: 



/ao 



'A 



/ao 



0(c 



/a{i + ^(Mac^) - d A c(tx)] - [d A (*2) - d A (h)f) + 

+ ^ ((v A • dc)(*a) - (v A • dc)(<i) + Y u ^b)dt) } + 0{Ar 2 , c" 3 ). 



(C20) 



Note that the vectors Rac (given in Eq. Ra and Rb, (given in Eq. ([791 ) are measured simultaneously at 

the time of signal reception. The result of Eq. (|C20I) is what we need in order to derive the DOWR and DWORR 
observables of the GRAIL mission given by Eqs. (|78p and With the X-band navigation on GRAIL with frequency 
~ 8.4 MHz, the result is accurate to ~0.5 mHz, which is sufficient for the purposes of navigating the GRAIL twins 
around the Moon. 

The expression for frequency transformation, Eq. (|C20|) . relates the frequency of the signal transmitted by the 
GRAIL spacecraft as reported at the BCRS to that measured by the DSN station. This result is complete up to the 
1/c 2 order, sufficient for GRAIL. 



